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ABSTRACT 

The work described in this dissertation can be divided into two parts. The first 
part is an investigation of the transient behavior and stability property of a phase 
conjugate resonator (PCR) below threshold. The second part is an experimental and 
theoretical study of the PCR’s spatiotemporal dynamics above threshold. 

The time-dependent coupled wave equations for four-wave mixing (FWM) in a 
photorefractive crystal, with two distinct interaction regions caused by feedback from 
an ordinary mirror, was used to model the transient dynamics of a PCR below thresh- 
old. Analytical expressions of the steady state cavity’s fields for the case of nonde- 
pleted pumps and an absorption free medium were derived and used to determine the 
self-oscillation conditions. The solutions, through simple frequency domain transfor- 
mation techniques, were used to define the PCR’s transfer function and analyse its 
stability. 

Taking into account pump depletion and medium absorption, the transient buildup 
and decay times of the cavity s fields as well as the specularly reflected and phase 
conjugate reflected intensities were numerically calculated as functions of a number 
of system parameters such as the coupling parameter and the pump and probe ratios. 

General trends were unveiled and discussed in view of the possible use of the PCR in 
image storage or processing architectures. Experimental results for the buildup and 
decay times confirmed qualitatively the predicted behavior. 

Experiments were carried out above threshold to study the spatiotemporal dynam- 
ics of the PCR as a function of the Bragg detuning achieved by misaligning one of the 
two pump beams and of the degree of transverse confinement controlled by varying 
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the resonator’s Fresnel number. The temporal aspects of the beam s complexity were 
studied by local intensity time series, power spectra, and reconstructed pseudo phase 
portraits. An irregular time series was identified as deterministic chaos having a cor- 
relation dimension of 5.2 and a Kolmogorov entropy of 0.16 s . Experimental results 
varying both control parameters revealed the presence of two distinct frequencies in 
the power spectra in some regions of the parameter space. The existence of optical 
vortices in the wavefront were identified by interferometry. The spatial complexity of 
the beam was studied in terms of the spatial distributions of optical vortices, their 
trajectories and their relationship to the beam’s spatial coherence. 

The transverse dynamics and the spatiotemporal instabilities were also described 
by modeling the three dimensional coupled wave equations in photorefractive FWM, 
using a truncated modal expansion approach. Numerical solutions of the model re- 
vealed the presence and motion of optical vortices in the wavefront. Simulations 
using the Bragg detuning as a control parameter showed that optical vortices appear 
as long as a large enough modulation of the wavefront is brought about by local gain 
enhancement resulting from the small misalignment of one of the pump beams from 
the Bragg angle. Power spectra of simulated time series also contain two distinct 
frequencies for a certain range of the off-Bragg parameter. Maps of the spatial cor- 
relation index were found to be templates of the corresponding maps of the vortices’ 
trajectories, indicating that the dynamics may be defect-mediated. 
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Chapter 1 


Introduction 


1.1 Background 

Optical feedback has been widely used in image processing since the late 1970s. 
Its applications include image restoration and recovery, pattern recognition, and the 
implementation of optical associative memory. Historically, optical feedback in ana- 
log processors was first proposed to extend their range of optical operations and to 
increase their performance and accuracy. Systems with passive coherent feedback 
have been demonstrated in spatial filtering, computing the solutions of differential 
equations, image restoration, and the implementation of nonlinear operation [1], Yet, 
the performance of these passive feedback systems is limited by the cavity losses and 
by the accumulation of phase errors in each round trip. 

Theorectical work on optical feedback suggested that fast analog optical devices 
with feedback can implement slowly converging iterative algorithms [2, 3], Itera- 
tive algorithms are often needed in problems such as image restoration and recovery 
as well as in recognition and associative recall. Efficient optical implementations of 
these algorithms, however, require active feedback with a regenerative feedback loop 
in which the signal energy and its phase after each iteration can be restored. In the 
feedback loop, to amplify the light, two type of gain devices have been experimen- 



CHAPTER 1. INTRODUCTION 

tally demonstrated in all-optical systems: dye amplification [4, 5] and photorefractive 
amplification [6, 7]. 

Amplified phase conjugation via four-wave mixing (FWM) with two antiparallel 
pump beams in a photorefractive medium can not only compensate for the cavity 
losses, but also cancel the aberrations and phase distortions acquired in each round 
trip. One advantage of photorefractive nonlinearities over other nonlinear mechanisms 
is that efficient beam interactions and large gain coefficients can now be achieved 
with a few tens of milliwatts of laser power in materials such as bismuth silicon oxide 
Bii 2 Si0 2 o- Low-power cw phase conjugate resonator (PCR) oscillations were recently 
demonstrated using a barium titanate BaTi0 3 crystal as a phase conjugate mirror 
(PCM) [8]. 

Since the successful demonstration of the photorefractive oscillators, a variety of 
all-optical architectures was proposed for the realization of optical neural comput- 
ers [9, 10, 11]. Most of the systems demonstrated so far are dealing with associative 
recall in which a partial information input can address the holographically stored 
images. The image having the highest correlation peak is selected and reconstructed. 
The PCMs in the systems are used to provide optical feedback, thresholding and gain, 
the three essential building blocks of a neural network system [12, 13, 14]. 

About a quarter of a century ago, an optical phenomenon now known as the 
photorefractive effect was discovered in lithium niobate LiNb0 3 and called at that 
time the “optical damage”. The undesirable optical damage in nonlinear and electro- 
optical crystals turned out to be an advantage when using these materials as holo- 
graphic optical memories. However, the interest in holographic memory in these 
media quickly declined because the stored information was being erased during the 
read-out process or lost subsequently after the pumps are turned off. For a while, as 
the basic understanding of the origin of the effect and its potential use in eliminating 
wavefront distortion became known, photorefractive wave mixing and phase conju- 
gation has grown rapidly into a sub-discipline in the field of nonlinear optics. In the 
middle 1970s, the science of the photorefractive effect became mature. More knowl- 
edge of the electro-optic media and their characteristic properties was available. The 
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previously unwanted degradation of optical memories was recognized as exactly what 
is needed for real time holography. Up until a few years ago, most of the research 
in that field was concerned with the steady-state behavior of wave mixing in pho- 
torefractive media and its applications [15, 16]. Recently, however, there has been an 
increasing interest in the dynamics of these processes. As more practical devices are 
built exploiting photorefractive phase conjugation, we would like to know how these 
systems respond to the temporal variations of any of their parameters. Knowledge 
of the dynamical behaviors of such systems makes it possible to design devices and 
predict their performance. Also of interest, even for steady-state applications, is the 
stability of various wave-mixing schemes. For a device based on these principles, one 
would like to predict if the system will be stable for the ranges of parameter values 
over which the desired operations must be performed. 

The dynamics of wave mixing in photorefractive media were first considered in 
early theorectical works of the transient behavior of two-wave mixing (TWM) [17]. 
They showed, for example, how the system approaches the steady state when only 
one of the input beams is present and the other is abruptly turned on. Later work 
examined the stability of photorefractive FWM processes [18, 19]. In addition, many 
groups have recently reported observations of temporal instabilities and chaos in 
photorefractive devices such as self-pumped PCMs [20, 21, 22], ring phase conjuga- 
tors [23, 24], mutually pumped phase conjugators [25, 26], and PCRs [27, 28]. Several 
models have been put forward to explain the dynamics of optical phase conjugation 
found in photorefractive crystals [29, 30, 31, 32] and devices [33, 34, 35]. Different 
mechanisms have been suggested to explain the temporal photorefractive instabilities. 
In the model of reference [22], for example, the chaotic behavior of a photorefractive 
PCM in a FWM geometry is shown to be caused by the existence of multiple inter- 
action regions in the crystal. Instabilities also arise in PCMs having only a single 
grating and a single interaction region [31]. Here, the factor leading to the chaotic 
behavior is the presence of an external electric field causing a shift in the optical fre- 
quency of the phase conjugate wave (PCW). A more recent model [32] even predicts 
a chaotic behavior in the standard FWM geometry with a single interaction region 
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and no external electric field if two gratings, transmission and reflection, of similar 
strength take part in the process. In ref. [35], two coupled ring cavities provide a 
competition between nonlinear gain and loss which results in nonsteady oscillating 
beams. It is thought that the competition for energy between different cavities re- 
sembles the competition for energy between different channels inside the self-pumped 
PCM [21] when it becomes unstable. 

Thus far, all these models are based on the well-known one dimensional coupled 
wave theory [36]. Therefore, they are incapable of describing the dynamics in the 
transverse plane observed in experiments such as those reported in refs. [27, 28]. 
Besides, effects such as strong light scattering (fanning) in the crystals, or various 
spatial effects of diffractions in resonators, all of which are found in real situations, 
result in non-uniform gain of the beam profiles. These kinds of considerations will 
require an extension of the one dimensional framework to a two or three dimensional 
framework in the coupled-wave theory. Indeed, the geometry of the self-pumped 
PCM has been considered in a three dimensional theory [37]. Furthermore, finite 
gaussian beams have recently been used to study the beam profile deformation in 
photorefractive TWM and FWM in a two dimensional framework [38]. These studies 
have been dealing with the steady-state aspect of the transverse beam profiles in 
photorefractive media. Yet, the transverse dynamics of nonlinear optical system is 
becoming an increasingly important topic of research [39]. 


1.2 Dissertation outline 

In this thesis, the study of the nonlinear dynamics in a PCR has been divided into 
two parts. In the first part, we will investigate the transient behavior and stability 
property of a PCR below threshold. In the second part, we will study experimentally 
and theoretically the PCR’s spatiotemporal dynamics above threshold. 

In chapter 2, we will present the model and the steady state analysis of a PCR 
below threshold based on time-dependent coupled wave equations for FWM in a 
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photorefractive crystal with two distinct interaction regions caused by feedback from 
a slightly tilted ordinary mirror. The steady state equations for the cavity’s field, 
taking into account nondepleted pumps and absorption free medium, will first be 
solved. The analytic solutions are used to discuss the threshold conditions. We also 
find the PCR s transfer function by applying simple frequency domain transformation 
techniques to the steady state equations. The transfer function is used to analyse the 
stability of the PCR. 

In chapter 3, the study will be extended to the transient regime of the PCR. We 
will first apply a numerical scheme for integrating the coupled wave equations using 
an adiabatic elimination process. We calculate the cavity’s buildup and decay times 
as well as the specularly reflected and phase conjugate reflected intensities as func- 
tions of a number of system parameters. Numerical results indicate that the cavity 
buildup and decay times can be tailored by varying several of the system parameters 
to assess the possible use of the PCR as an image processing element. Experimental 
measurements are shown to qualitatively confirm the numerical simulations. 

In chapter 4, we will give a brief introduction to some of the terminlogy and 
methods of nonlinear and chaotic dynamics, which are applied to the analysis of 
experimental time series data and to the beam’s spatial complexity in the dynamics 
of PCR above threshold. 

In chapter 5, we will present experimental results for the dynamics of PCR above 
threshold. We study by experiments some aspects of the spatiotemporal behaviors of 
the PCR as a function of the resonator’s Fresnel number by varying one of the two 
apertures inside the cavity. Temporal aspects of the beam’s complexity will be studied 
analysing the local intensity time series, power spectra, and embedded phase space 
portraits. A range of parameters is found for which the time series exhibit chaotic 
variations. For the spatial aspects of the beam’s complexity, we will give evidence for 
the presence of wavefront dislocations or vortices in the optical field by using video 
recording and interferometry. A rich variety of defect movements are recorded and 
analysed. The spatial complexity is characterized by the number of optical vortices 
and their variety of movements as well as their relationship with the beam’s spatial 
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coherence. 

In chapter 6, we will present and test a model based on a truncated modal expan 
sion of the cavity modes, to investigate the spatiotemporal dynamics of the PCR. We 
find the optical vortices in the optical field’s wavefront in the numerical solutions. Nu- 
merical calculations are performed using the Bragg detuning, achieved by misaligning 
one of two pump beams as a control parameter. Numerical results indicate that the 
loss of temporal coherence and the onset of temporal chaos in the local intensity fluc- 
tuations are correlated to the loss of spatial confinement of the vortices’ trajectories 
and the spatial coherence. The validity of the model is verified by comparison with 
experimental data. 

In chapter 7, a conclusion and a summary for this work are given. Some possible 
future investigations are suggested. 


1.3 The photorefractive effect, wave mixing and real time 
holography 

In the remainder of this chapter, we begin a brief introduction to the photorefrac- 
tive effect, wave mixing and phase conjugation, and then go on to the transients of 
the PCR and the steady state analysis in the next chapter. 

The photorefractive effect is a phenomenon in which the local index of refraction 
is changed by a spatially varying light intensity. The effect has been observed in 
crystals such as lithium niobate LiNbOs, barium titanate BaTiC >3 and strontium 
barium niobate SBN as well as in cubic crystals of the sillenite family such as bismuth 

silicon oxide BSO. 

Suppose that we illuminate a photorefractive crystal with two coherent plane elec- 
tromagnetic waves as in fig. 1.1a. The waves will interfere with each other and pro- 
duce a spatially varying (sinusoidal) light intensity pattern inside the crystal. Charge 
carriers are optically excited into the conduction band in the bright regions of the 
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Figure 1.1: Pictorial illustration of the formation of a phase hologram in electro- 
optic crystals via the photorefractive effect, (a) Two coherent electromagnetic plane 
waves write a sinusodial grating inside the crystal, (b) Spatial variations of the light 
intensity I(x ), charge carrier distribution p(x), space-charge field E(x) and refractive 
index n(x). After Hall et al. [15] 
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interference pattern and move by drift (due to any internal and/or externally applied 
electric fields) and diffusion into the dark regions, where they are trapped. Thus, 
this sinusoidally varying space charge distribution creates a static, spatially periodic 
electric field E(x) inside the crystal as shown in fig. 1.1b, according to Poisson’s 
equation: 


V • E(x) = (4w/c)p(x). (1.1) 

Here p(x) is the charge density, e is the dc dielectric constant. The electric field E(x) 
alters the crystal lattice and locally modulates the crystal’s refractive index n(x) via 
the linear electro-optic effect, generating a phase diffraction grating. (See fig. 1.1b). 
The amplitude A n(x) of the change in index of refraction is [40] 


An(x) = —(n 3 r eJ j/2)E(x). 


( 1 . 2 ) 


Here r cSf is the effective Pockels coefficient of the crystal. (Pockels coefficients are lin- 
ear combinations of elements of the crystal’s third rank electro-optic tensor Rijk [40].) 
Crystals having large Pockels coefficients such as BaTi0 3 and SBN can give rise to a 
significant photorefractive nonlinearity. 


In a seminal series of papers, N. V. Kukhtarev et al. [41, 42, 43] showed that this 
refractive index grating interacts with the light waves and couples them together, so 


that their intensities 1 \ and I 4 and their phases </> 1 and <j ) 4 change as the light waves 
propagate through the crystal: 


dR 

dz 
<U 4 
dz 

d<t>\ 

dz 

d(f>4 

dz 


= — Re( 7 0 ) ~r~7 ~ a h, 


Re(7o)- 


I\ + I4 

hh 


I1 + I4 

= 


a/ 4 , 


(1.3) 


Here a is the absorption coefficient of the crystal, and 7 0 is the coupling strength 
(see eq. B.54), which is primarily determined by the crystal’s Pockels coefficients, 
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dielectric constant and interacting geometry. In the absence of any external fields 
and if the charge migration in the crystal is dominated by diffusion, the resulting 
index grating will be shifted by 90° with respect to the interference pattern, as shown 
in fig. 1.1b, and the coupling strength 7 0 will be purely real. 

It is the occurrence of the phase shift between the index grating and the interfer- 
ence pattern that makes possible an asymmetric transfer of energy from one beam 
to the other. This two-beam coupling phenomenon is known as two wave mixing. 
According to equation 1.3, if the real part of the coupling strength, Re( 7 G ), is greater 
than zero (it depends on the interaction geometry and the sign of the charge of the 
charge carrier), the intensity / 4 will grow with propagation while I\ will drop. Also, 
for purely real 7 0 , the phases of both beams remain unchanged as they travel through 
the crystal. 

Suppose now that a third weak beam (beam 2) is incident to the crystal, in a 
direction opposite to that of beam 1 as shown in fig. L2. This beam will be diffracted 
by the phase grating into the direction opposite to that of beam 4. It can be shown 
that this diffracted beam (beam 3) is the phase conjugate of beam 4. To see what 
this means physically, let the beam 4 amplitude be written as 

E 4 (r,t) = ~ [A 4 {r)e i( - t? ~^ + A\{r) e -'^-^] , (1.4) 

The phase conjugate wave (PCW) is obtained by taking the complex conjugate of 
the spatial part of the wave in eq. 1.4: 

E pc (r, t)= l - + A 4 (r)e~^- tr -^] , (1.5) 

which can be rewritten as 

E P c(r, t) = + A 4 (r)e’<‘ V+ "'>], (1.6) 

which is the time reversed version of the wave in eq. 1.4. 

This generation of a time-reversal copy of the probe wave is called optical phase 
conjugation. A photorefractive crystal employed to produce such a phase conjugate 
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Figure 1.2: Schematic diagram of the four- wave mixing. Beam 1 and beam 4 write a 
phase grating inside the crystal. Beam 2, counterpropagating to beam 1, is incident 
on the crystal and diffracted into the direction opposite to that of beam 4. The 
diffracted beam, beam 3, is a phase conjugate replica of beam 4. 


10 


CHAPTER 1. INTRODUCTION 


wave back to its source is called a “phase conjugate mirror”. The process we just 
described involves four beams interacting and is referred to as “four- wave mixing”. 

H. Kogelnik, more than 25 years ago, demonstrated using conventional holography 
that a distorted image could be undistorted by first producing its phase conjugate 
replica and then sending it backward through the original distorter [44], In conven- 
tional holography, this involves a two-step procedure in which a photographic plate is 
first used to store the interference pattern between the incoming distorted image and 
a reference wave. This hologram is developed and repositioned exactly in its original 
place. The hologram is then read out by a wave propagating in a direction opposite 
to that of the original reference wave and deflects the reading wave to produce the 
PCW. 

Now, with the process of photorefractive FWM, the photorefractive crystal re- 
places the photographic plate to store the hologram which is being written by an 
image-bearing beam and a coherent reference beam. At the same time, the hologram 
can be read out by a reading beam incident at the Bragg angle inside the crystal. 
The simultaneous write and read processes eliminate the developing and fixing steps. 
So the FWM gives us a new method of storing dynamic holograms in real time image 
processing, which is often named “real-time holography”. 
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Chapter 2 


The PCR model: steady state 
analysis 


2.1 Introduction 

Phase conjugate resonators (PCRs) are resonators in which at least one of the 
mirrors is a phase conjugate mirror (PCM). The PCM, which acts both as a phase 
conjugator and as a gain element, confers to the resonator some unique properties 
(e.g. phase healing, transverse mode degeneracy, etc.) [45]. A number of applications 
of the PCR to image processing have been demonstrated or proposed [7, 46, 9, 11]. 
A class of these applications makes use of the regenerative feedback of the PCR to 
synthesize unusual transfer functions, perform nonlinear operations, or implement 
iterative algorithms [7, 46]. Another type of problem which has motivated this study 
more directly involves the implementation of discrete time algorithms for distributed 
parameter control systems [47]. A PCR could be used as a real-time holographic 
element for temporal image storage and delay operations. In all these, the PCR is 
below threshold (for self-oscillation) and functions as a regenerative amplifier. 

In these applications, the way in which the PCR responds to time- varying signals is 
important for understanding the operation of the device. Therefore a theory is needed 
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to describe the transient effects in the PCR (such a theory is described in 2.2). Phase 
conjugate Fabry-Perot cavities have recently been studied both in the steady state 
case [48] and in the limit of a PCM with a fast response such as a Kerr material [49]. 
Fast response means that the time scale associated with the buildup of the polarization 
field in the medium is much shorter than other time scales in the problem. Thus, 
there is no difficulty in coupling the time dynamics of the material with the time 
dynamics of the optical field in a self-consistent manner. However, this is not the 
case in sluggish photorefractive materials such as in BaTi0 3 because the time scale of 
the grating formation is not shorter than other time scales in the problem. Equations 
for the transient behavior of the PCR with a sluggish medium must, therefore, include 
the material equations. 

In section 2.3, we first solve the steady state equations for the PCR, for the con- 
ditions of nondepleted pumps and an absorption free medium. An analytic solution 
can be obtained and used to describe the threshold conditions of the PCR. In the 
quasi-steady state limit for a slow material, the steady state solution can also be used 
to derive a transfer function for the PCR (section 2.4). The transfer function is then 
used to analyse the transient response if there are some system operative ranges that 
could lead to instability of the device. 


2.2 Model of the PCR below threshold 

In this section, we develop a model for a phase conjugate Fabry-Perot resonator 
with a PCM consisting of a sluggish photorefractive material. The medium is modeled 
by the standard Kukhtarev theory. Equations for the fields in the cavity and in the 
medium, under the plane wave theory, are derived. 

2.2.1 The Cavity 

The geometry of the PCR is sketched in Figure 2.1. It consists of a conventional 
(dielectric) mirror and an externally pumped PCM . The input beam A IN may make 
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Figure 2.1: Schematic diagram of a phase conjugate resonator using an externally 
pumped phase conjugate mirror (PCM) and a conventional mirror (M). r m , t m , r' m , t' m 
are the reflection and transmission coefficients of the conventional mirror, fi i, fi 2 
are the PCM reflectivities at the two interaction regions. The angle of incidence 
of the input beam and the distance between the two interaction region are greatly 
exagerated 
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a small angle with the mirror normal so that in general two distinct interaction 
regions exist in the PCM. The presence of two distinct gratings can easily be verified 
experimentally. 

A PCR formed by two PCMs has been investigated by others [50] and proposed 
for the holographic storage of images [51]. The only difference between the cavity 
of fig. 2.1 and the two-crystal cavity, as shown in figure 2.2, is that, in the two- 
crystal set up, the pumps’ intensities can be varied independently. In a single-crystal 
cavity where the two interaction regions share the same pumps, the steady state PCM 
reflectivities (i 4 and \i 2 are equal, although they generally grow at different rates. 

Within the framework of a plane wave theory, the fields’ amplitudes in the cavity 


(-L < z < 0), satisfy the self-consistency equations 

A 4 i(z,t) = t m A IN (z,t) + r' m A 32 (z,i - ^ e'^ 2( - L+z \ (2.1) 

A 12 (z,t) = r m A 3 1 (z,t - ^ e *7 2 ( L +*) } (2.2) 

A 3 i(z, t ) = Hi ft + 0 A* 4l (z, t + , (2.3) 

A 32 (z, t) = /i 2 (i + A 42 ^ 2 , t + — ^ , (2-4) 


where L is the optical length of the cavity. 

The PCM is assumed to be a sluggish photorefractive medium such as BaTi 03 . 
When the medium response time is much longer than the cavity round trip time, all 
the temporal delays in equations 2.1 to 2.4 can be discarded. 

The specularly reflected output and the phase conjugate output are respectively 
Ar{-L, t) = r m A IN {-L , t) + CA 32 (0, t)e‘c L , (2.5) 

and 

A c (-L,t) = t' m A 31 (0,t)e^ L . (2.6) 


15 



CHAPTER 2. THE PCR MODEL: STEADY STATE ANALYSIS 



Figure 2.2: Two-crystal cavity comparable to the cavity with two interaction regions 
of figure 2.1. 
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2.2.2 The Four-wave mixing (FWM) equations 

The standard model used to describe photorefractive wave mixing is based on the 
coupled wave formalism originally developed in the late 1960’s by Kogelnik [36], and 
on the model due to Kukhtarev et al. [42] which describes how the index grating forms 
in the photorefractive medium. We first outline the coupled wave equations for the 
photorefractive FWM (the coupled wave equations are derived in Appendix A) and 
then combine them with Kukhtarev’s equations to give a set of coupled partial differ- 
ential equations which completely describe both the temporal and spatial variations 
of the fields in the FWM process. 

The beams in the PCM are treated as plane waves with slowly varying envelopes 
(SVE approximation); i.e., the temporal and spatial rates of change of the fields and 
material polarization amplitudes are slow compared to the light frequency and coarse 
compared to its wavelength. The pump and input beams have all the same frequency 
and the same polarization state. 

With these assumptions, the wave equation in the medium ( 0 < z < d ) leads to 
the following standard coupled equations (see Appendix A) (the subscript j = 1 or 2 
refers to either one of the two interaction regions): 

( - ajAifat), (2.7) 

f d d \ 

~ y^Qf- J t) = ~ Qj( Z i 0^3 j( z i 0 + a jA 2j( z > (2-8) 

f d d \ 

~~ vdt J ^ = Qj( z, > t) + OL 3 Azj{z, t), (2.9) 

(J; + A 4 j( 2 >t) = Qj(z,t) A v(z,t)- 0 ' j A^(z,t), (2.10) 

where cxj is the linear amplitude absorption coefficient, Q 3 is the complex amplitude 
of the material grating and v is the speed of light in the medium. The angular 
dependence of the spatial derivatives has been neglected, which is consistent with the 
assumption of a small angle between the pump beams and the cavity axis. 
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If the medium response time is much longer than the transit time through the 
PCM, as is assumed here, the temporal derivatives can be dropped from the LHS of 
equations 2.7- 2.10. The equations then describe a quasi steady state situation in 
which the fields’ amplitudes follow adiabatically the medium changes. 

Futhermore, it was assumed that only one index grating is operative in both 
regions. For example, with a geometry utilizing the extraordinary polarization to 
take advantage of the large electro-optic coefficient of BaTiOs, the dominant grating 
is the one formed by the interference of Ai with A\ and A 2 with A3. Its wave number 
is ki — k 4 = k 3 — k 2 , where k 3 is the wave vector of beam j. 

The dynamics of the grating formation can be derived from the Kukhtarev equa- 
tions and has the following form (see Appendix B where the details of the Kukhtarev 
theory and its underlying assumptions are summarized): 

Tj— Qj( z ,t) = Qoj( z ) ~ Qj( z it)i j = (2-11) 


where Tj is a possibly complex time constant and Q 0 j is the steady state grating. 
These two quantities can be evaluated from the parameters in the Kukhtarev model 
as shown in section 2.2.3. From Appendix B, the steady state grating has the form 


Qoj{z) 7oj 


AuWA^jz) + AlMAzjjz) 


7 = 1,2, 


( 2 . 12 ) 


where I OJ is the local average intensity in region j : I 0 j(z,t) — 52 \Akj(z,t)\ 2 - The 
second factor in the RHS of equation 2.12 is the modulation depth of the interference 
fringes of beam A\ with A 4 and beam Ai with A3. 7oj I s coupling coefficient and 
has the form 


loj — 


iijn 3 0 

4c 


E’sc , 


' effuse 


(2.13) 


where is an effective Pockels coefficient which depends on the geometry and the 
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beams’ polarization, and E s sc is the complex amplitude of the space-charge field mod- 
ulation (see Appendix B). The effective Pockels coefficient r e// for barium titanate 
(point group 4mm) can be written [7] as 

r <// = ^r 13 cos/?, (2.14) 


for ordinary polarized beams, and as 


r e/f — 


n o r i 3 (cos 20 cos 2/9) -f 4n 2 rc 2 r 42 sin 2 /9 + rc^r^cos 2/2 + cos 2#)] cos /9, 


(2.15) 

for extraordinary polarized beams, where n 0 and n e are the respective ordinary and 
extraordinary indices of refraction at wavelength A = 27rc/o;, are the elements of 
the linear electro-optic tensor (Pockels coefficients) [40], and 0 and /3 are the angles 
inside the crystal defined as in Figure 2.3. 


If the small angular difference for the beams interacting in regions 1 and 2 is 
neglected, 7 0 and a have the same values in both interaction regions. 


2.2.3 The material’s constants 

According to the Kukhtarev band transport model, a mathematical description of 
the grating formation process is given by the equations 


dn_dN* 1 , 

dt dt q 1 

(continuity), 

(2.16) 

dt ~ ^ + - Nq) ~ IrkNc, 

(rate equation), 

(2.17) 

J = q\mE — k B TpVn + ple c 

(current equation), 

(2.18) 

V-E=l(n + N A -N+) 

(Poisson’s law), 

(2.19) 


where 
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C-axis 


t 



Figure 2.3: Geometrical configuration of two- wave mixing in a photorefractive mate- 
rial. K is the wave vector of the induced grating; c is the optical axis of the crystal, 
20 is the angle between two beams inside the crystal; /? is the angle of K with the 
c-axis. 
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n(r,t) 

N d 

NMrJ) 

N a 

n d -n+ 

j 

s 

I 

1r 

p 

E 

k B 

T 

P 

9 


is the free charge carrier (single species) 
number density, 

is the total number density of the dopant (single species) 
(Fe for example), 

is the number density of the ionized donors (traps) 

(Fe 3+ for example), 

is the number density of the acceptors, 

is the number density of neutral (filled) 

donors (Fe 2+ for example), 

is the current density, 

is the photoionization cross section, 

is the intensity of the optical wave, 

is the recombination rate coefficient, 

is the mobility of the charge carriers, 

is the total local electric field, 

is Boltzmann’s constant, 

is the temperature, 

is the photovoltaic coefficient, 

is the charge of the charge carriers, and 

is the static dielectric constant. 


After linearization in the grating modulation index and in the steady state limit, 
equations (2.16)-(2.19) can be solved and the steady state space-charge field is found 
to be (see Appendix B): 


where E 0 is the external 
material: 


E ? 


—iE„ 


( E 0 + Ep) + iE B 


E 0 + i(E q + E D ) J’ (2 ' 20) 

applied field, E p , E D and E q are characteristic fields of the 


E P 

E d 


pL 

qpn 0 

k B T\K\ 

9 


(photovoltaic field), 
(diffusion field), 


( 2 . 21 ) 
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E q = 


qNe 

clKI 


(saturation space charge field). 


Here I a is the zero order of the light intensity, N e is the effective trap number density, 
K is the grating vector and p is the photovoltaic coefficient [52]. 


With no external field (E a = 0) and negligible photovoltaic field (E p ~ 0), the 
steady state space-charge field is simply given by 


_ -lEgEp 

sc E V + E D 


( 2 . 22 ) 


In this case of a diffusion dominated process, the index grating in equation B.52 
leads the light interference fringes by a 7 t/ 2 phase shift and the coupling constant 
7 „ in equation 2.13 becomes real. The model also predicts a complex response time 
constant of the form (see Appendix B) 


To = T e + i , 




(2.23) 


where the material response time is 


1 

(>+5) 

2 

1 + 

fla 

\TE 

y 

0 + ^) 0 + 

ia) 

TD ) 

+( 

s) 2 (?f) 


and the material frequency is 


(2.24) 


U 


9 


i te)^- 1 ) 
M' + s) 2 +(“) 2 ' 


(2.25) 


The various time scales in equations 2.24 and 2.25 are defined by 


6 (dielectric relaxation time), 

q\in Q 

(2.26) 

t , = (n\K\E„)-' 

(drift time), 

(2.27) 

= (c|k|£d)-' 

(diffusion time), 

(2.28) 
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t r — (TflA^r 1 (recombination time), 

Tj = ( SI 0 ) 1 (inverse of photoexcitation rate), 
and the free carrier average density is 


n Q = 


SI 0 (N d -N a ) 

1rN a 


(2.29) 

(2.30) 


(2.31) 


All the material constants involved in these relationships are of great interest to the 
material physicist whose objective is to make materials with optimal properties [53], 
It is apparent that different applications require different material properties. For 
instance, in order to demonstrate high speed, low write laser power and long-lived 
memory, or large gain, it is necessary to use several different types of materials. It has 
been observed that BSO is a fast photorefractive material but of low efficiency due 
to a small electro-optic coefficient while BaTi 03 has by far the largest electro-optic 
coefficient leading to high diffraction efficiency but with a slow response time. The 
research for improving photorefractive materials has a goal to overcome the trade off 
between speed and sensitivity as well as to enlarge the materials’ performance over 
a wider operative range of device exploiting these materials, which is essential to the 
new optical computing technologies. 

For photorefractive materials, one important attribute of the the grating response 
time is its dependence on the total intensity incident on the crystal. The one-species 
model of the Kukhtarev theory predicts an inversely linear dependence (r oc I~ l ). In 
BaTi0 3 , this dependence is experimentally found to be sublinear with an exponent 
0.5 < x < 1. This behavior can be accounted for by assuming that two different 
species of traps take part in the photorefractive process [54]. For our study, it is 
however sufficient to use a normalized response time of the form 


( N a \ 1 A 

r - ~ T - = ^ kfrnrj t. = v (2 - 32) 

which is valid when the recombination time t r is short compared to the other time 
scales of the model. This happens to be the case for BaTi0 3 . In the calculations of 
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the next chapter, we used a normalized time scale r Q = I D 1 • This normalized time 
scale can be related to the real time scale for a particular material by evaluating the 
constant A. Typical numbers for BaTi0 3 are given in Ref. [53]: 


7 h = 

N a ~ 

fi ~ 

S(Nd — Na) ~ 

9 = 


5 x 10 -8 cm 3 /s, 
2 x 10 16 cm -3 , 
0.5cm 2 /V ■ s, 
a ~ 0.3cm _1 , 
1.6 x 10- 19 C, 


e 0 = 8.85 x 10 14 Fcm *. 


(2.33) 


The relative dielectric constant is given by 

e R = C|| sin 2 0 + cos 2 /3, (2.34) 

with e|| = 4300, e x = 168 and 0 is the angle between the grating wave vector and the 
c-axis of the crystal (Figure 2.3). For the experiment we had 0 ~ 10°, which gives 
t R ~ 300. With these numbers, we find A = 1.1 x 10 18 cm 2 . 

In equation 2.32, I 0 is expressed in number of photons per cm 2 per sec. In the 
calculations of the next chapter, the unit of intensity was chosen to be 1 mW/mm 2 . 
At a wavelength A = 0.514/im (photon energy = 3.87 x 10- 19 J), this corresponds to 
a photon flux of 2.58 x 10 17 s _1 cm -2 . The normalized time scale used is thus given by 


with 


r 0 (Normalized) = [/ 0 (in mW/mm )] 


r(in sec) 

A' 


1.1 x 10 18 
2.58 x 10 17 


4.26(s mW/mm 2 ). 


(2.35) 

(2.36) 


2.3 Steady-state analysis 


The main reason for seeking an analytic steady state solution to the PCR s equa- 
tions is that in the quasi steady state case (slow medium such as BaTi0 3 ) such a 
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solution can be used (as shown in section 2.4) to define a transfer function from 
which the stability of the solution can be analyzed. Furthermore, the steady state 
solutions determine the threshold for self-oscillations of the PCR. 

2.3.1 Steady-state equations and their solutions 

The steady state equations are obtained simply by setting all the time derivatives 
equal to zero. At steady state, an analytic solution can easily be found in the case 


of nondepleted pumps and negligible absorption, i.e. A\{z) — Ai(0) = A\, A 2 (z) = 
A 2 (d) = A 2 , and a\ t 2 = 0- Equations ( 2.7- 2 . 12 ) then reduce to 

A/l 3) ( 2 ) = j + A-A:Jz)\ A 2t (2.37) 

2-A- ij (z) = ^[A 1 Al J (z) + A- 2 A 3j (z)\A\ j = 1,2, (2.38) 

where 

h~\Atf+\A 2 \ 2 . (2.39) 

The boundary conditions are, from equations ( 2 . 1 - 2.4), 

^4i(0) = 0) + r , rn As2{0)e Xc2lj ', (2.40) 

2l«(0) = 4A 3 ,(0)e‘““, (2.41) 

A3i{d) = 2432 (d) = 0. (2.42) 


To solve the steady state equations 2.37 and 2.38, one may proceed as follows. 
Multiply equation 2.37 by p = A\j A x and we have 

— (pA 3 j(z)) = lA-A- A* 4 j (z) + pA 3j (z)^ , p = (2.43) 

j-A^z) = To (l - fa(z) + P A 3 j (z)} , (2.44) 
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Adding equations 2.43 and 2.44 and integrating, gives 

pA 3l (z) 4 A' tj (z) = Mai(O) + A^.(0)]e w . (2.45) 

Substituting this equation into equation 2.43 and 2.44 and integrating, gives, respec- 


tively 

\ A 1 2 



= LA [p4 3j (0) + A^(0)| (e y °* - 1) + ^(0), 

(2.46) 

and 




AjjW = [p4 3j (0) 4 ^(0)] (e y ° z -1)4 A,,(0). 

(2.47) 


Using the boundary conditions 2.42 in equation 2.47, we find 


^(0) = 

Substituting this in equation 2.40 and 2.41, we find 


R v 


h 


£-7 od _ | 

A 4 i( 0) = * m A/jv(0) + r' m p * , — A4 2 (0)e'- 2L 


and 


^42(0) = r' m p 


e~^ d + R p 

e -io d _ i 


e~' li + R , 


-A 4 i(0)e 


i — 2L 


Equations 2.49 and 2.50 are combined to give 

tmAiwiO) 


with 


and 


P = \r m \ 2 R P 


A 4 i(0) = 

-lod _ ^ |2 


e-i° d + R p 


R P = 


A\ 


\-p ’ 

(feedback parameter), 
(pump ratio). 


(2.48) 


(2.49) 

(2.50) 


(2.51) 


(2.52) 


(2.53) 


Using this result in equations 2.46 and 2.47 and making use of the relationship 
between A 42 (0) and A 41 (0) derived in equations 2.49 and 2.50, we finally find the 
steady state fields in the crystal (0 < z < d): 
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A 4l( Z ) — 


e 7o( *~ rf) + / t m A IN (0) ' 
e~' 1 ° d + R p l 1-/3 


(2.54) 


Ki( z ) = r mPe 


-^2L (g 7o(2 ~ d) +^)(e~ 7od -ir ( LnAlN^Y 

\e~^ d + R p \ 2 l 1-/3 


(2.55) 


^31 (*) — P 


e 7o( *~ d) -l / <m^w(0) y 
e -7o<i + R p \ 1-/3 / 


(2.56) 


^ 32 ( 2 ) = r'^Rpe 


.i* 2 L {e^ z - d) - l)(e~' y ° d - 1)* (t m A IN (0)' 


|e -7od + R p \ 2 


1-/? 


(2.57) 


The cavity fields and the outputs can then be found by using equations ( 2.54 
2.57) in equations ( 2.1- 2.6). 


2.3.2 Threshold for self-oscillation 

The PCM steady state reflectivities are defined, from equations 2.3 and 2.4, as 

4jj(0) 


Pi 


^(o y 


j — 1 , 2 . 


(2.58) 


Using equations ( 2.54- 2.57), we find 


. e~ 7od - 1 
Pi = Pi = P = P — 


e -io d -j- R P ’ 


(2.59) 


and the feedback parameter takes the form 


/3 = \r m p | 2 - 


(2.60) 


If the two interaction regions share the same pumps, they have equal steady state 
reflectivity. In a diffusion dominated medium with a stationary grating, the coupling 
parameter 7 0 is real. The phase of the PCM reflectivity is thus simply the sum of 
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Pump ratio R p 


Figure 2.4: Steady state PCM reflectivity as a function of the pump ratio for various 
values of the coupling parameter 
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the phases of the pump beams: Arg[jz] = <^(0) + <fa(d), <f>i( 0) = Arg[Ai(0)], fa(d) = 
Arg[A 2 (d)]. 

Figure 2.4 shows the magnitude of the PCM reflectivity as a function of the pump 
ratio for various values of the coupling constant 7 Q d. For the pump ratio R p <C e~ lod , 
\fi\ increases as yjRp. For R p > e~ lod , |//| drops as 1 / yjfy, almost independently of the 
coupling constant. Also, the reflectivities all converge to zero as the power of pump 2 
becomes much greater than that of pump 1. This is because the index grating written 
by pump 1 and the input is erased by the large power of pump 2. Equation 2.60 shows 
that the feedback parameter is proportional to \fi\ 2 . Care must be taken, however, 
in interpreting this result. It is clear that the steady state solution of equations 
( 2.54- 2.57) is valid only below threshold, i.e. /? < 1. The condition /9 = 1, which 
corresponds to a feedback system in which the gain exactly compensates the losses, 
defines the threshold for self-oscillations of the PCR (i.e. no external input needed 
to establish the cavity fields). Figure 2.5 shows /? calculated from equation 2.60 with 
a mirror reflectivity |r m | 2 = 0.95. It is seen that for small coupling (7 a d ~ 1.8), the 
resonator is stable for the entire range of pump ratios. For larger coupling constant, 
there is a range of pump ratios around R p ~ 0.1 which would have /3 > 1. This is 
of course not physical. As one approaches threshold, the cavity fields increase, the 
undepleted pump assumption ceases to be valid, saturation takes place and /? remains 
equal to unity even if the coupling constant is increased. 


2.4 Transfer function of the PCR and its stability analysis 

In this section, we develop a method to derive a transfer function of the PCR and 
we use it for the analysis of the device’s frequency response and possible instabilities. 

2.4.1 Frequency domain transformation of the PCR’s equations 

To derive the transfer function, we begin with the time-dependent FWM equations 
( 2.7- 2.11) in section 2.2.2. In the undepleted pumps and negligible absorption 
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CCL 



Figure 2.5: Feedback parameter as a function of the pump ratio for various values 
of the coupling parameter. The steady state solution is valid only below threshold: 
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approximation, they read: 

Gk + Am) = QA z ,t) A 2, 


' d_ _ _d_' 

dz vdt 


A 4 j( zA ) = Qj( z H) A A 


(~ 4\ A \ A lj{ zA ) + A 2 A 3j( z , 0 

T Qj \ z ’ 0 — 7 o j 


(2.61) 

(2.62) 

(2.63) 


We define a 3 j(z, il), fi 4 j(z,Sl) and qj(z, ft) as the complex Fourier transforms (or 
double sided Laplace transforms) of A 3 j(z,t), A* 4j (z,t) and Qj{z,t ), respectively, i.e. 

r+oo 

(2.64) 

(2.65) 


/"TOO 

F(t)e-' nt dt, 

-oo 


1 /* + co+tc 

f (<) = 5- / my m dn, 

Z7T «/ — oo+tc 


so that a 4 (z,17) = a 4 (z,— f)). Equations ( 2.61- 2.63) become 


d 

.ft 


+ l — 

dz 

V 

d 

ft 

dz 

V 


^4 j(-2,17) — Qj Li) , 

(1 + iQr)qj(z,Q) = y[Aia A j{z ,D) + A* 2 a 3 j(z, Q)]. 

*0 

Thus, the complex amplitudes of a 3j and a 4j satisfy 

( d Q\ i m 1° A \a 4j {z,-D) + A;a 3i (z,fi) 

(al + 'v) °* (z ' n) = mn; Z 


(2.66) 

(2.67) 

( 2 . 68 ) 

(2.69) 


'0 ^ * / m _ 7 0 + ^a 3j (z,fi) 

dz i v r 4AZ ’ U) -\+iDt z ! - 


(2.70) 


For a slow medium, the transit time through the PCM is much shorter than 
the medium response time (r » d/v) and the terms ifl/v can thus be neglected. 
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Equations 2.69 and 2.70 are then identical to the differential equations 2.37 and 2.38 
used to derive the steady state solutions if the substitution 


Jo 

1 + iDr ’ 


(2.71) 


is made. With a slow medium (r > cavity round trip time), the boundary conditions 
for 03 and a 4 also have the same form as equations ( 2.40- 2.42). The solutions of 
equations 2.69 and 2.70 are thus identical to the solutions of equations 2.37 and 2.38 
respectively, with the substitution of equation 2.71. 


2.4.2 Transient response and its stability 


From the results of the previous section, we can now define a transfer function for 
the phase conjugate cavity beam, 


with 


11(D) 


<*31 (0; fl) 
«/jv( 0 , D) 

( t m p \ 

\i -m) 


e~ ld - 1 

e-^ + R,' 


(2.72) 


m 



e~ ld - 1 2 
e~^ d + R p ’ 


(2.73) 


where p — A* 2 j A\ and R p = \ AzlA\\ 2 . 

From the definition of the complex Fourier transform given in equations 2.64 
and 2.65, the system is unstable when the transfer function has poles with negative 
imaginary parts. Figures ( 2.6- 2.8) show contour plots (on a log scale) and 3D plots 
(on a linear scale) of the transfer function for three different values of the coupling 
parameter. The shaded areas outline the region in which /? > 1 and in which the 
definition of a transfer function becomes meaningless. The poles are located along 
the boundary j3 = \. 

For a weak coupling parameter, all the poles are above the real axis and the sys- 
tem is stable (Figure 2.6). From equation 2.73 and the parameter values given in 
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Figure 2.6: (a)Contour plot in the complex plane (log scale) and (b)3D plot (linear 
scale) of the transfer function of the PCR with a sluggish medium, in the limit of 
negligible absorption and no pump depletion. The relevant parameter values are 
7 0 d = 10, R p = 0.1, R m = 0.95. The shaded region outlines the domain f3 > 1. 
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Figure 2.7: Same as Figure 2.6 with 7 0 d = 2.0. The cavity has just reached threshold. 
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- 2-1012 



Figure 2.8: Same as Figure 2.6 and 2.7 with 7 Q d = 3.0 showing additional unstable 
conjugate pairs of poles. 


35 




CHAPTER 2. THE PCR MODEL: STEADY STATE ANALYSIS 

the caption of Figure 2.6, it can be calculated that the intersection of the bound- 
ary /3 = 1 with the imaginary axis crosses to the lower half plane for 7 a d = 1.87 
(Figure 2.7). This is, of course, precisely the coupling parameter needed to reach 
threshold for these parameter values, as can be verified in Figure 2.5. At this point, 
the resonator self-oscillates and has a stable steady state output. As the coupling 
parameter increases further, additional complex conjugate poles may cross the real 
axis and become unstable. (Figure 2.8). 

2.5 Conclusion 

We have presented a model and the steady state analysis of the PCR in the 
framework of the plane wave representation. A steady state analysis of the FWM 
process of the PCR is performed and analytic solutions are derived in the case of 
nondepleted pump, an absorption free medium and real coupling constant. Stability 
analysis of the solutions is also discussed with the use of the transfer function in 
fjre frequency domain. We find that dynamic instabilities start to come in when the 
resonator self-oscillates, i.e. when the feedback parameter f3 = 1. 
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Chapter 3 

Transient behavior of a PCR: 
numerical study 


The theory of the photorefractive phase conjugate Fabry-Perot cavity formulated 
in chapter 2 was used to find steady state solutions. The analytic results and the 
analysis of the stability under quasi steady-state conditions have been presented in 
that chapter. For practical applications, it is also crucial to understand the transient 
dynamics of the device. In this chapter, the transient regime of a PCR under the 
plane wave assumption is studied numerically. In particular we are concerned with 
the cavity’s response time as a function of the system’s parameters. In section 3.1 we 
review the time- dependent equations for the waves in a PCR and describe a numerical 
method for their solution. In section 3.2 we present the numerical results as functions 
of the coupling parameter, the pump ratio, the probe ratio, the dielectric mirror’s 
reflectivity, and the absorption constant of the medium. The general trends are 
discussed in view of the possible use of the system in image processing. In section 3.3 
we describe an experimental setup used to study the transient dynamics of a PCR 
and discuss the experimental results. A summary of the study is given in section 3.4. 
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3.1 A numerical approach 


In this section we begin with the time-dependent equations derived in chapter 2, 
including the pump depletion and medium absorption. We then apply an adiabatic 
elimination method to integrate the partial differential equations separately in the 
space and time domains to obtain the final steady state of the cavity fields. 

3.1.1 The start of the transients 

The time-dependent equations derived in section 2.2.2 for the plane wave model 
of four wave mixing in a slow photorefractive medium are (figure 2.1). 

—Aij(z,t) = — Qj( z i t)A-4j{ z i I) — a jAij(z, t), 

oz 

>(*,<) = Qi(z,t)A,j(z,t) + ajA 3 j(z,t), 

OZ 
OZ 

and 

r j t Q,(z,t) + *) = 0 + A- 2 ,(z,t)A 3 ,(z, ()] 

The boundary conditions for the cavity fields are: 


(0,f) = t m AlN + r m^32(CM) e * c2L > 

(3.6) 

A 42 {0,t) = r' m A 31 (0,t)e^ 2L , 

(3.7) 

A3i(d, f) — ^32 (d> — 0, 

(3.8) 

-An (0, i ) = -4i2(0,t) — Aio, 

(3.9) 

^2l(d,0 = ^22 {d,t) — ^2o* 

(3.10) 


(3.1) 

(3.2) 

(3.3) 

(3.4) 

(3.5) 
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For the initial state of the PCR, we assume that only the two pump beams Ai and 
A 2 are present. We further assume that there is no scattered noise so that even if the 
cavity is above threshold (gain > loss), self-oscillations do not start spontaneously 
without an input beam. Since the interaction between the two pumps is neglected, 
their initial spatial distributions in the crystal assume a simple exponential depen- 
dence of the form e ±az . 

The signal beam is then abruptly turned on to a constant value Ain- Since no 
interaction takes place yet, it also assumes an exponential distribution: 


A 4l (z,0) — t m Ai N e az , (3.11) 

^3i(*,0) = A 42 (z,0) = A 32 (z, 0) = 0. (3.12) 

After a first short lapse of time, a grating is being formed in region 1, driven by the 
interference between A 41 and A x only. The scattering of A 2 by this grating produces 
a phase conjugate beam A 3 i- After reflection on the mirror, this beam produces a 
probe beam A 42 in the region 2. The three beams A 42 ,A\ and A 2 in region 2 vary 
exponentially in space since there exists no interaction yet, and their interference 
triggers the formation of a grating in exactly the same way as in region 1. After a 
second lapse of time, this grating generates a phase conjugate beam A 32 which adds 
up coherently to the input Ain to define the new probe beam A 44 in region 1. The 
process goes on with the two gratings driving each other’s growth and building up 
the cavity fields until a steady state is reached. 


3.1.2 Adiabatic elimination algorithm 

For a slow medium such as BaTi0 3 , the time derivatives in equations ( 3. 1-3.4) 
have been eliminated since the propagation delays in the cavity and in the medium are 
short compared with the time needed for grating formation. Under this quasi steady- 
state condition, one can regard the field amplitudes inside the crystal as following 
the slow grating changes adiabatically. This method was proposed by Heaton et al. 
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in [17] for the integration of the time-dependent photorefractive equations separately 
in the space and time domain and is known as the adiabatic elimination method. The 
numerical recipe includes the following steps [55]. At time tw, Q modulates 
the field amplitudes inside the crystal. The spatial distributions of the fields are 
calculated by solving the ordinary differential equations ( 3. 1-3.4) with the boundary 
conditions (3.6-3.10), by using a fourth order Runge-Kutta method. The boundary 
conditions allow one to separate these equations into two independent sets. First, 
equations 3.2 and 3.3 are solved simultaneously with the boundary values A 3 j(d, t N ) = 
0; A 2j (d,t N ) = A 2j ( 0) = A 2o . The solution of this first set of equations is then used 
to find the boundary values A 4 j (0, t^). These, together with the boundary values 
Aij(0,tjv) = Aij(0) = Aio, are used t° solve equations 3.1 and 3.4. All the fields 
inside the crystal at are then known. 

The medium then responds to these fields according to equation 3.5 to produce 
a new grating Qj(z,tN+ 1 )- A second order Runge-Kutta method is used to integrate 
this temporal equation, giving 


Qj(z,tn-n) S3 Qj[z,tn) + (A t K /T){i a mj(z,t N + At N /2) - [Qj(z,t N ) + (l/2)/\ ( ' v) ] } , 

(3.13) 


where 


/C W) = (At lV /T)['f 0 m ] (z,t N ) - <2,(z, /fv)], 


(3.14) 


and 


A\j(z, f yv) "f A 2 j{Z) In) 


(3.15) 


with tjv+i = tN + A tN. The new grating is then used to find the new spatial distribu- 
tion of the fields at time t N +i which is again calculated from equations 3.1 and 3.4. 
The process is repeated for a fixed number of steps or until a steady state is reached. 

For the calculation to be valid, the time step A tjy must be small enough compared 
to the time constant T(f;v) of the grating. Since this constant is a function of the 
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intensity incident on the crystal and therefore changes in time, the temporal step 
was not taken as a constant but rather each step was chosen equal to 1/10 the 
instantaneous time constant r(tjv) of the grating. 

Since, according to the standard Kukhtarev model, the grating time constant is 
inversely proportional to the total intensity I 0 incident on the crystal, we chose a 
relative time scale of the form r = I~ x for the calculations. This time scale can 
be related to the actual time scale T — AI~ X by evaluating the constant A for a 
particular material [53]. Typical values of parameters for BaTiOa give A ~ 4.3 s 
x mW/mm 2 (see section 2.2.3). The time appearing in the graphs discussed in the 
next section is normalized to the medium response time with a total incident I a of 1 
mW/mm 2 (e.g., ~ 4.3s for BaTiOa). 

3.2 Numerical analysis 

In this section, the numerical procedure just discussed is used to find the transient 
behaviors of the cavity fields as functions of the system’s parameters. Particular 
attention is given to the system response time and the PCR’s reflectivity because of 
their importance in the design of practical devices. 

3.2.1 Response curves for the cavity fields 

The system’s parameters used in this and the following sections are the coupling 
parameter 7 0 d, the pump ratio R p = I 2 /I 1 , the probe ratio R s — // ; v / ( I\ + h), the 
input mirror reflectivity R m = |r m | 2 , and the absorption parameter ad. Figures ( 3.1- 
3.5) show the results of sample calculations of the temporal evolution of the cavity 
field intensity / 3 i( 0 , t/r 0 ). The other fields experience similar growth. The normalized 
intensity of the phase conjugate cavity field intensity / 3 i/(/i + I 2 ) is plotted as a 
function of time, for a square pulse input. The length of the pulse was chosen long 
enough to allow the PCR to reach a steady state. 
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Figure 3.1: Normalized cavity field intensity / 3 i/(/i + h) showing the response of 
the PCR to a square pulse input for various values of the coupling parameter and 
R a = 0.02, R p = 0.3, R m = 0.95, ad = 0.15. Time is normalized to the response time 
t q of the medium with I Q = lmW/mm 2 . The cavity is above threshold for 7 Q d > 2.5. 
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Figure 3.1 shows the response for various values of the coupling parameter. The 
other parameters are R s = 0 . 02 , R p ~ 0.3, R m = 0.95 and ad - 0.15. Larger values 
of the coupling parameter lead to larger cavity fields, as expected, to faster growth 
rates and to shorter buildup times. For 7 0 d ~ 2 . 2 , the cavity is above threshold and a 
self-oscillating field remains in the cavity even after the input has been cut off. This 
is consistent with the result of figure 2.5 which indicates that at R p ~ 3, the cavity 
reaches threshold for 7 Q d ~ 2 without absorption. 

Figure 3.2 shows the response for various pump ratios. In varying R p the pump 
intensity /1 was maintained constant, and I 2 was varied. The other parameters are 
7 0 d — 2 , R s = 0 . 02 , R m — 0.95 and ad — 0.15. With these parameter values, the 
cavity remains below threshold for all values of R p , as shown in figure 2.5. It is seen 
that at first, the steady state cavity field intensity increases with increasing pump 
ratios and then drops. This is consistent with the result of figure 2.4 which shows 
that the PCM reflectivity also reaches a maximum. The buildup and decay rates 
follow similar trends but the total buildup time does not vary much with the pump 
ratio. 

Figure 3.3 shows the PCR response for various probe ratios. The other parameters 
are R p — 0.5, 7 0 d = 2.5, = 0.95, and ad = 0.15. The steady state cavity intensity 

increases monotonically with the probe ratio, the buildup rate increases and the 
buildup time becomes shorter. The decay time, however, is fairly constant. In this 
example, the cavity was slightly above threshold. As the input is cut off, the cavity 
field evolves toward a final self-oscillating steady state intensity in a time which does 
not depend much on the initial input intensity. This behavior could be found useful 
in an optical processor, as it is seen that seeding the PCR with inputs of various 
intensities changes the buildup time but does not affect the threshold condition. 

Figure 3.4 shows the response for different mirror reflectivies and for 7 0 d = 2, 
R p — 0.5, R s = 0.1. In this case, lowering the input mirror reflectivity increases the 
steady-state PCM reflectivity product Consequently, the increased PCM gain 

overcomes the additional losses of the input mirror and leads to a larger steady state 
cavity field. The buildup and decay times slowly increase with smaller input mirror 
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Figure 3.2: Normalized cavity field intensity / 3 i/(/i + h) showing the response of 
the PCR to a square pulse input for various values of the pump ratio and 7 a d = 2.0, 
R a = 0.02, R m = 0.95, ad = 0.15. 
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Figure 3.3: Normalized cavity field intensity / 3 i/(/i + / 2 ) showing the response of 
the PCR to a square pulse input for various values of the probe ratio and 7 Q d = 2 . 5 , 
R p = 0.5, R m = 0.95, ad = 0.15. 
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Figure 3.4: Normalized cavity field intensity I 3 i/(h + h) showing the response of the 
PCR to a square pulse input for various values of the input mirror’s reflectivity and 
lo d = 2.0, R p = 0.5, R s = 0.1, ad = 0.15. 
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Figure 3.5: Normalized cavity field intensity / 3 i/(/i + / 2 ) showing the response of the 
PCR to a square pulse input for various values of the medium’s absorption parameter 
and -f 0 d = 2.0, R v = 0.5, R a = 0.02, R m = 0.95. 
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reflectivity, since this brings the cavity closer to threshold. 

Figure 3.5 shows the response for different values of the absorption parameter. 
Increasing the absorption losses reduces the cavity field amplitude, everything else 
being equal. Here, 7 0 d = 2, R p = 0.5, R s = 0.02, and R m = 0.95. The buildup and 
decay rates also drop with increasing absorption. 

Most of the response curves shown in figures ( 3. 1-3.5) exhibit smooth and mono- 
tonic growth to steady state and monotonic decays. With large coupling constant and 
especially if the resonator is far above threshold, overshooting and damped oscillations 
can also be seen. 


3.2.2 Cavity response time 

The cavity buildup and decay times were measured from a large number of plots 
similar to those shown in figures ( 3. 1-3. 5). These times were then plotted as functions 
of the pump ratio and the probe ratio, which ratios appear to be the two most 
important parameters in determing the system’s behavior. 

Since some of the buildup curves exhibit overshooting for large 7 Q d, the buildup 
time was defined as the time it takes for the cavity intensity to grow from 10% to 90% 
of its maximum (which is not necessarily the steady state in case of overshooting). 
The decay time was defined as the time it takes for the cavity intensity to drop from 
90% to 10% of its steady state value. 

Figures 3.6 and 3.7 show the buildup and the decay times, respectively, as functions 
of the pump ratio for (a) R m = 0.8, (b) R m = 0.95 and R s = 0.1. For weak 
coupling, the buildup time increases with increasing pump ratios for small R p , reaches 
a maximum and drops for larger R p . This can be correlated with the behavior of the 
PCM reflectivity and the cavity field intensity which follow similar trends with varying 
R p , as seen in figures 2.4 and 3.2. For stronger coupling, there is a range of values of 
R p within which the cavity is above threshold. In this region, the initial gain of the 
cavity is large (larger than the losses) and the transient buildup is fast, sometimes 
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Figure 3.6: a Buildup time versus pump ratio for various values of the coupling 
parameter and R m = 0.8, R s = 0.1, ad = 0.15. For 7o d > 2.5, the cavity is above 
threshold for a range of values of the pump ratio around 10 -1 . b Same as figure 3.6a 
except that R m = 0.95. 
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exhibiting overshooting. This explains the central dip in the curves (figure 3.6) for 
7 od > 2. 

In figure 3.7, the decay time is plotted as a function of the pump ratio. The 
divergence of the decay time as threshold is approached is clearly shown. The pump 
ratios at which divergence occurs correspond to the ratios at which the feedback 
parameter /? reaches unity in figure 2.5. 

Figures 3.8 and 3.9 show respectively the buildup time and the decay time as 
functions of the probe ratio for various values of the pump ratio and for 7 Q d = 2.5, 
R m = 0.95, ad = 0.15. For these values of parameters, the resonator is above thresh- 
old for pump ratios ranging from about 0.02 to 0.5. Far below threshold (R p ~ 1), 
figure 3.8 shows that the buildup time does not change much for a wide range of probe 
ratios. As one approaches threshold ( R p « 0.01 and 0.55), the buildup time varies 
significantly with the probe ratio especially for small probe beam intensity. Figure 3.9 
shows a similar behavior for the decay time. Below threshold, the decay time does 
not change with the probe ratio while close to threshold, it varies significantly with 
that parameter. This may bear some important consequences in applications making 
use of the cavity for image storage or processing. If a large dynamic range input is 
imaged in the PCM, the weaker parts of the input will take considerably more time 
to buildup and to decay than the more intense parts, giving rise to contrast enhance- 
ment during buildup and contrast reduction during decay. The phenomenon could 
also be used to produce momentary time varying spatial filtering if the input spatial 
spectrum rather than the input itself is projected in the PCM. 

To simulate the buildup of self-oscillation in a cavity above threshold, we calcu- 
lated the cavity response to a very weak input pulse of short duration. In the physical 
system, this type of seed is provided by scattered noise. The result is shown in fig 
ure 3.10. In figure 3.10a, the cavity is below threshold. A very weak phase conjugate 
beam appears and remains in the cavity as long as the seed is present. It then rapidly 
decays when the seed is removed. In figure 3.10b, the cavity is above threshold. The 
cavity field remains in the cavity even after the seed is removed and it grows very 
slowly (note the different intensity and time scales of figure 3.10a and b). After a 
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Figure 3.7: a Decay time versus pump ratio for various values of the coupling pa- 
rameter showing the range of values of pump ratio for which the cavity is above 
threshold. The parameter values are the same as for figure 3.6a (R m = 0.8). b Same 
as figure 3.7a but for the same parameter values as for figure 3.6b (R m = 0.95). 
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Figure 3.8: Buildup time versus probe ratio for various values of the pump ratio and 
lo d = 2.5, ad = 0.15, R m = 0.95. 
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Figure 3.9: Decay time versus probe ratio for various values of the pump ratio and 
the same parameter values as for figure 3.8. 
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relatively long time, the cavity field builds up more rapidly to a final steady state 
value. 

3.2.3 Specular and conjugate reflections 

We now use the numerical method described in the preceding sections to investi- 
gate the temporal response of the two reflected outputs, I R = \A R \ 2 and I c = |Ac| 2 . 
Figures 3.11 and 3.12 show the two outputs normalized to the total pump intensity 
/j + / 2 , as functions of time for a normalized input R s — 5 x 10 -3 . The results of fig- 
ure 3.11 correspond to a pump ratio of 0.5 and a modest coupling parameter 'y 0 d 2.5. 
Three different values of the input mirror reflectance R m were used (figures 3.11a-c). 
It is seen that with a smaller input mirror reflectance, the total PCM steady state re- 
flectivity product |/xi// 2 | increases, approaching unity. Consequently, the steady state 
specular reflection almost drops to zero, leaving only the phase conjugate reflected 
wave counter propagating to the input beam. 

With a higher coupling parameter (figure 3.12), the PCM steady state reflectivity 
exceeds unity. In this case (figure 3.12a) the specular reflection is seen to drop to 
zero and then to settle at a finite steady state value which is small compared to 
the conjugate reflection. The situation is that of a cavity which is not too far from 
threshold. Similar temporal conjugate reflection was reported in [49] for a PCM 
having an instantaneous response and a fixed reflectivity. Decreasing the input mirror 
reflectance again increases the PCM reflectivity, bringing the cavity above threshold 
(figure 3.12b). In this case, the specular reflection drops to zero and then grows to a 
large steady state value. 

3.3 Experimental results 

An experimental setup was built to study the temporal response of the PCR as a 
function of the pump ratio, the probe ratio and other parameters. Measurements of 
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Figure 3.10: Response of the PCR to a short seed pulse of intensity corresponding 
to R s = 10 -24 . Normalized cavity field I 31 versus normalized time for R p = 0.1, 
Rm = 0.95, oid = 0.15. a Cavity below threshold with / y 0 d — 2.0; b Cavity above 
threshold with 7 Q d = 2.2. (Note the different scales). 
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Figure 3 . 11 : Time evolution of the normalized intensity of the specular Ir/{I\ + h) 
and conjugate Ic j (-/i T I 2 ) reflection with a constant normalized input of R s 0.005, 
R p = 0.5, 7 0 d — 3 0. a, R m = 0.95 resulting in a steady state PCM reflectivity 
product \'n lf i 2 \ = 0.915; b, R m = 0 . 8 , |/ii/z 2 | = 0.921; c, R m = 0 . 6 , = 0.939. 
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Figure 3.12: Same as figure 3.11 with larger coupling parameter 7 Q d = 6.0. a R m = 
0.95 resulting in I/11/Z2I — 1-008; b R m = 0.6, resulting in l/^i^l = 1-444 which brings 
the cavity above the threshold. 
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the cavity’s response times are described and compared, qualitatively, with theoretical 
predictions in this section. 

3.3.1 Experimental setup 

A diagram of the experimental setup is shown in figure 3.13. The PCM is an 
externally pumped single crystal of BaTi0 3 . The laser is a single mode Ar + , operating 
at A = 514nm. The two pump beams and the input beam are each equipped with 
a halfwave plate/polarizer assembly used as a variable attenuator. All polarizations 
are extraordinary (in the plane of the figure). The pump and the probe beams made 
angles of about 49° and 15°, respectively, with the normal to the crystal surface. The 
input mirror (dielectric) had a reflectivity of 95%. It made a small angle with the 
input beam in order to produce two distinct interaction regions in the crystal. 

The phase conjugate cavity beam was extracted from the cavity by a pellicle beam 
splitter PBS having about 5% reflectance. The lens Le reimaged the interaction 
regions in the plane of a slit which was positioned to select the cavity beam /31 
and isolate it from the other beam I 3 2 . The beam intensity I 3i was measured by a 
photodetector, recorded and stored in a PC. A shutter was used to turn the input on 
and off and produce a square pulse, as in the simulation. 

3.3.2 Measurement of the cavity response time 

Two sets of data were collected - the first for as high a coupling constant as we 
could achieve with this particular geometry and the other for a lower coupling con- 
stant, spoiled by rotating the crystal by a few degrees. In both cases, the response to 
the input pulse was recorded and the buildup time and the decay time were estimated 
from the curve. A typical experimental record of the cavity field intensity / 3 1 is shown 
in figure 3.14. 

The results are shown in figures ( 3.15-3.18). The error bars in these figures 
are qualitative and are intended only to give an indication of the level of confidence 
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Figure 3.13: Diagram of the experimental setup. L: Ar + laser (514nm), A/2 -P: 
halfwave plate/polarizer assembly, Sh: shutter, M: 95% dielectric mirror, PBS: pellicle 
beam splitter, Le: imaging lens, S: slit, D: photodetector. 
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Figure 3.14: Typical cavity field response to a square pulse input from which the 
buildup and decay times were measured. 
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that was attributed to the eyeball estimation of the characteristic times from the 
response curves. Some data were noiser than others and thus were assigned larger 
error bars. The lines are visual aids. Above threshold, the PCR response often 
exhibits overshooting and sometimes damped oscillations. This made it more difficult 
to measure the characteristic times, and larger errors are expected as threshold is 
approached. 

There were many other possible sources of errors more difficult to quantify. For 
example, close to threshold the system is extremely sensitive to any external pertur- 
bation. In this region, small mechanical or acoustical perturbations can give rise to 
large differences in the PCR transient response. 

Figure 3.15 shows the buildup time as a function of the pump ratio. Pump l x was 
maintained constant at 0.5mW with a beam diameter of about 1 mm inside the crystal, 
and pump I 2 was varied. The probe ratio was R s = 0.02. The general behavior is 
qualitatively similar to that of figure 3.6 obtained from the simulation. At low pump 
ratios, the buildup time increases with R p and with the coupling parameter. It reaches 
a maximum for 0.1 < R p < 0.2 and then drops. At high pump ratios, the buildup 
time is nearly independent of the coupling parameter. The two lower data points near 
the center of the high coupling curve are dubious. They could represent the drop of 
the buildup time predicted in figure 3.6 above threshold or their lower value could be 
due to some external perturbations giving rise to large errors near threshold. 

Figure 3.16 shows the decay time as a function of the pump ratio. This also is 
qualitatively similar to the result of the simulation shown in figure 3.7. By comparing 
the range of pump ratios within which the experimental PCR is above threshold with 
the curves of figure 3.7b, which corresponds to the same values of parameters, one 
can estimate that the coupling parameter achieved in the experiment was of the order 
of 7 0 d ~ 2 . 

Figures 3.17 and 3.18 show the buildup and decay times, respectively, as functions 
of the probe ratio for 7 Q d ~ 2, R p ~ 1 . The buildup time decreases monotonically 
with R s in a manner quite similiar to that shown in figure 3.8. The decay time 
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Figure 3.15: Experimental buildup time versus pump ratio for two values of the 
coupling parameter and R m = 0.95, R s = 0.02. The error bars only indicate the level 
of confidence to be given to different data points. The solid lines are for visual aid. 
(The circles are data points for a low coupling parameter; the squares are data points 
for a high coupling parameter) 
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Figure 3.16: Experimental decay time versus pump ratio for two values of the coupling 
parameter and R m = 0.95, R s = 0.02. (The circles are data points for a low coupling 
parameters; the squares are data points for a high coupling parameter) 
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behavior is also qualitatively similar to that shown in figure 3.9. 


3.4 Conclusion 

In this chapter we have presented a numerical scheme for integrating the pho- 
torefractive equations separately in the space and time domains using the adiabatic 
elimination process. With this numerical scheme, the transient dynamics of a PCR 
was described in the limit of a plane wave model and a sluggish photorefractive 
medium, but taking into account both the pump depletion and medium’s absorption. 
The dependence of the cavity buildup and decay times on parameters such as the 
pump ratio, the probe ratio, the coupling parameter, the input mirror’s reflectivity 
and the medium absorption constant were discussed. The numerical results show that 
the parameters affecting the threshold most directly are the coupling parameter and 
the pump ratio. We find that the higher the coupling parameter, the wider the range 
of pump ratio values for which the PCR is above threshold. Results also indicate 
that the cavity buildup and decay times can be tailored by varying several of the 
system parameters. This may be of importance to image processing applications. For 
example, the PCR would, in principle, perform contrast enhancement during buildup 
and contrast reduction during decay. This same feature could be used to implement 
a time varying spatial filter if the projected field in the PCM is the spatial spectrum 
of an image rather than the input image itself. 

The numerical simulations also reveal the unusual reflection property of a phase 
conjugate Fabry-Perot. For example, we find that when the steady state value of the 
phase conjugate reflectivity approaches unity, the specular reflection is completely 
extinguished, leaving only the phase conjugate reflection output. 

An experimental setup was built to study the transients in a PCR and the exper- 
imental results are in general qualitative agreement with the results of the numerical 
simulations. In real situations, the material’s parameter values may vary considerably 
between crystals. In our case, the parameter values for the BaTi0 3 crystal are cho- 
sen from the literature [53]. Futhermore, it seems likely that quantitative agreement 
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would require some refinement of the simplified Kukhtarev model by including, at 
least, multiple species of traps taking part in the photorefractive process. 
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Temporal and spatial instabilities 


So far we have demonstrated how the PCR below threshold responds to the tem- 
poral change of an input. In the remaining part of this thesis we will study another 
mode of operation of the PCR: emphasis will be given to the characterization of the 
spatial structure of the cavity beams above threshold. In this chapter, we first give 
a brief introduction to some of the terminology and methods of chaotic dynamics, 
which will be applied to the study of the spatiotemporal dynamics of the PCR above 
threshold in the following two chapters. The first part of this introduction will be 
devoted to chaos in nonlinear dynamics. This includes the discussion of qualita- 
tive methods in nonlinear dynamics such as the uses of power spectra, phase space 
portraits and the notion of attractors. Regarding the quantitative methods for the 
characterization of chaotic dynamics, we give a brief introduction to the notion of 
dimensions and to the calculation of the correlation integral function as well as the 
Kolmogorov entropy. We also describe a method used for the reconstruction of a 
pseudo phase space portrait from an experimental time series. In the second part, we 
will briefly mention some recent studies of spatial pattern evolutions and instabilities, 
which is currently an active area of investigation [56]. The discussion will focus on 
wavefront dislocations and optical vortices because they are found to appear in the 
PCR’s beams and probably mediate the observed dynamics. 
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4.0.1 Nonlinear systems 

In general, the dynamics of physical systems can be described by means of dif- 
ferential equations. In cases of nonlinear systems, the equations which govern their 
behaviors will contain nonlinear terms. Linearization of the nonlinear terms in some 
regions of parameter space where their effects are small often provides simple solutions 
and a sufficiently clear picture of the system’s dynamics. A classic example of this 
is the simple pendulum in the small angle approximation. There are some nonlinear 
systems, however, where the nonlinear terms dominate in the physically interesting 
region of parameter space. In those cases, the linearization of the system’s equations 
results in the loss of significant knowledge about the system’s behavior. The variety 
of responses of nonlinear systems is often found to be much richer than that of linear 
systems, as illustrated in figure 4.1 for two systems driven by a sinusoidally vary- 
ing input. The output of the linear system must also be sinusoidal, with the same 
frequency as the input. The output in the nonlinear system can be periodic (not 
necessary sinusoidal); or periodic with subharmonics with frequencies u/n, where n 
is an integer. Another possibility is the occurrence of quasi-periodic motion with an 
output of the form cicos(uqf) + C 2 cos(cu 2 t), where the ratio uq/uq is an irrational 
number if uq and u? 2 are incommensurate; or the output can be chaotic, fluctuating 
erratically in time. 

Now, let us look at one example of nonlinear system which is at the origin of 
the science of chaos. Chaos was first discovered by E. N. Lorenz [58] in the early 
1960s. Lorenz considered a rather typical problem of fluid dynamics which is called 
the convection instability or Rayleigh-Benard instability (his purpose was to analyze 
the motion of the atmosphere). The motion of the fluid is described by the Navier- 
Stokes equations which are nonlinear, partial differential equations. In order to cut 
down the complexity of the problem of solving these equations, Lorenz introduced 
a Fourier decomposition. He expanded the velocity and temperature fields of the 
fluid into spatial Fourier series. The Fourier coefficients were still time dependent 
variables. Lorenz truncated the infinite series to only three terms representing three 
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Periodic 

Subharmonic 

Quasi-Periodic 

Chaotic 


Figure 4.1: Schematic diagram showing the possible outputs of a linear and a non- 
linear system for a sinusoidal input. In a linear system, the output must be also 
sinusoidal. In a nonlinear system, the output can be either periodic, periodic with 
subharmonics, quasiperiodic or chaotic. (After [57]) 
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individual modes to save computing effort while still keeping the important physics 
of the problem. In this way he derived three coupled differential equations for three 
variables. The equations of the Lorenz model are usually written in the following 
form: 


dX (V Y\ 

-dT = a{Y - x) ’ 

(4.1) 

= -XZ + rX -Y, 
dt 

(4.2) 

dZ 

— = XY - bZ. 

dt 

(4.3) 


X represents the amplitude of the fundamental velocity mode; Y represents the fun- 
damental temperature mode and Z represents the temperature second harmonic, a is 
the Prandtl number, r = Rf R c is the normalized Rayleigh number relative to the criti- 
cal Rayleigh number R c at which the onset of convection occurs, and b = 47 r 2 / {'K 2 -\~k \ ), 
where k\ is a dimensionless wave number. The equations 4. 1-4.3 are of quite simple 
structure. They are ordinary differential equations and contain only two nonlinear 
terms in XZ and XY . To the great surprise of many physicists and mathematicians, 
Lorenz found that the numerical solutions for some parameter values (e.g. a — 10, 
b — 8/3, r = 28) can exhibit, in addition to stationary and periodic states, irregular 
behaviors. Figure 4.2 illustrates the rapid divergence, characteristic of chaotic behav- 
ior, of two solutions which start from initial conditions differing only by about 1 part 
in a thousand. Figure 4.3 shows the plots of X versus Y and X versus Z where the 
trajectories stay most of the time around certain zones although the trajectories are 
never closed. 


4.1 Characterization of chaotic dynamics 

4.1.1 Fourier spectra 

One method of examining the dynamics of a nonlinear system is by Fourier spectral 
analysis. Consider the temporal behavior of a function y(t). It is quasiperiodic if its 
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Figure 4.2: Graphs of chaotic fluctuations in time of the variable X of the Lorenz’s 
convective equations. Values of parameters are: o — 10, b — 8/3, and r = 28. This 
shows a chaotic solution very sensitive to initial conditions, i.e. two points which 
start out close can grow expontentially apart in time. 
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Figure 4.3: Trajectories of a chaotic solution of Lorenz’s convective equations pro- 
jected on the plane X — Y (above) and on the plane X — Z (below). The trajectories 
trace out certain regions of the planes although they are never closed. (Parameters’ 
values are the same as in figure 4.2.) 
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Fourier transform consists of sharp spikes, i.e. if 

= < 4 - 4 ) 

i= i 

The characteristic recurrence feature of quasiperiodicity may be stated formally as 
follows: for any e > 0 there exists a T(e) such that any interval of length T(e) of the 
real line contains at least one point t' such that \y(t) — y(t')\ < e for any t [59]. Given 
y(t), we can always find a t' such that y{t') is as close to y{t) as we wish, and there 
are an infinite number of such times if. Periodic functions are quasiperiodic, but of 
course quasiperiodicity does not imply periodicity. Thus the function 

y (t) = Ci cos u>it + c 2 cos u> 2 t (4-5) 

is quasiperiodic, but it is not periodic unless uq and u 2 are commensurate frequencies, 
i.e. unless ujiju} 2 is a rational number. 

Quasiperiodic motion can certainly look very complicated and irregular, but it 
cannot be truly chaotic in the sense of exponential sensitivity to initial conditions. In 
particular, the difference between two quasiperiodic trajectories is itself quasiperiodic, 
and so we cannot have the exponential separation of initially close trajectories that 
is the hallmark of chaos. 

Since quasiperiodicity implies order, it follows that chaos implies non-quasiperiodic 
motions. Thus chaotic motion does not have a purely discrete Fourier spectrum as 
in 4.4, but must have a broadband, continuous component in its spectrum. Fourier 
analysis is therefore a very useful tool in distinguishing regular from chaotic motion. 

Caution has to be taken here because a broadband spectrum can arise from 
stochastic behavior as well as from deterministic chaos. Additional analyses are 
therefore needed. For example, for a process governed by stochastic mechanisms, 
the spectral power at high frequencies follows a power law. In contrast, for a process 
governed by deterministic mechanisms, the power spectrum falls off exponentially [60]. 
Therefore, a linear behavior on the semilogarithmic plot of a power spectrum indi- 
cates that the broadband noise arises from a deterministic rather than a stochastic 
process. 
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4.1.2 Phase space portraits and attractors 

A dynamical system consists of two parts: the state which describes the essential 
information about the system, and the dynamics which is a rule that describes how the 
state evolves with time. The evolution can be visualized in phase space, an abstract 
space whose coordinates are the components of the state. Therefore, the phase space 
provides a powerful tool for inspecting the behavior of chaotic systems. The phase 
space represents the behavior in a geometric form. For example, a pendulum that 
moves with friction eventually comes to a rest, which in phase space means that the 
orbit approaches a stationary point called a fixed point. Since it attracts nearby 
orbits, it is known as an attractor. 

Some systems do not come to rest in the long term but cycle periodically through 
a sequence of states. An example is the pendulum clock (driven pendulum), in which 
energy lost to friction is supplied by an external force. The pendulum clock repeats 
the same motion over and over again. In phase space such a motion corresponds to 
a cycle, or periodic orbit. Such attractors are called limit cycles. 

A system may have several attractors. If that is the case, different initial conditions 
may evolve to different attractors. The set of points that evolve to an attractor is 
called its basin of attraction. The pendulum clock has two such basins. A very small 
displacement of the pendulum from its rest position results in a return to rest. With 
a larger displacement, however, the clock begins to tick as the pendulum executes a 
stable oscillation. 

Another type of attractor in phase space is a torus or quasiperiodic attractor, 
which resembles the surface of a doughnut. This shape describes motion made up of 
two independent oscillating frequencies such as in driven electrical oscillators. One 
frequency determines how fast the orbit circles the doughnut in the short direction 
and the other regulates how fast the orbit circles the long way around. If motion is 
represented by the combination of more than two oscillating frequencies, the attractor 
becomes a higher-dimensional torus. 

The important feature of a quasiperiodic attractor is that in spite of its complexity 


75 



CHAPTER 4. TEMPORAL AND SPATIAL INSTABILITIES 

it is predictable. Even though the orbit may never exactly repeat itself, i.e if two 
frequencies are incommensurate, the motion still remains regular. Orbits that start 
on the torus near one another remain near one another, and long term predictability 
is guaranteed. 

Since the first example of a chaotic or strange attractor discovered by Lorenz, now 
known as the Lorenz attractor, interest in looking for this special type of attractor 
in phase space has been growing rapidly. A chaotic attractor has a much more 
complicated structure than a predictable attractor such as a fixed point, a limit cycle 
or a torus. For a chaotic attractor, two orbits with nearby initial conditions diverge 
exponentially and then come back close together for only a short time. Stretching and 
folding operations which take place in phase space are the key to understanding the 
chaotic behavior on these attractors. Exponential divergence is a local feature because 
attractors have finite size so that two orbits on a chaotic attractor cannot diverge 
exponentially forever. As a result, the attractor must fold over onto itself. Although 
orbits diverge and follow increasingly different paths, they eventually must pass close 
to one another again. The process of stretching and folding happens repeatedly, 
creating folds within folds indefinitely. A chaotic attractor is, therefore, much like 
a fractal: an object that reveals more details but self similar as it is increasingly 
magnified. 


4.2 Dimensions 

To distinguish deterministic chaos from stochastic noise, other quantitative means 
are needed. A chaotic attractor can be quantitatively characterized by its geometric 
properties. A suitable quantity characterizing the attractor as a metric structure is 
its dimension. 

There are several alternative ways of defining the dimension of a set of points [61]. 
We shall give what we feel to be the most intuitive one here, called the fractal or 
capacity dimension D c [62]. Referring to figure 4.4, we consider a set of points 
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Me 



Figure 4.4: Illustration of the method of covering a 1- and 2-dimensional sets of points 
using small squares. This can lead to an intuitive definition of the capacity dimension 
(see the text). (After [57]). 
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uniformly distributed along a line as in the above half of the figure. We want to 
find out how many small squares of side length e are needed to completely cover the 
set of points. If the line is of length /, then this number will clearly be N(e) = l/e. 
Similarly, we can consider a two-dimensional set of points as shown in the bottom half 
of the figure, and find that the number of small squares needed to cover the area A 
is N(c) = A/e 2 . Extrapolating this procedure to sets of points in higher dimensions, 
we find the following scaling rule: 


l d 

N ^) = ^ 


(4.6) 


where d is the dimension of the set of points. From this we can define the capacity 
dimension Dc as 


Dr 


log N(e) 


log / + log(l/e) 

and in the limit of small e, the term containing l becomes negligible, i.e. 

logJV(e) 


Dr 


lim 


(4.7) 


(4.8) 


,C '-~6log(l/c)- 

If the dimension of a set of points as defined in equation 4.7 is non-integer, then the 
set is a fractal set. 


The structure of the strange attractor does not tell us anything about the fre- 
quency with which different points on the attractor are visited in the course of the 
flow. It is therefore useful to modify the definition of eq. 4.8 by covering the attractor 
with cells and weighting the jth cell by its visitation probability Pj. This gives the 
information dimension 

Dj = lim 1(e)/ log(l/e), (4.9) 

c — ►O 

where 

N(e) 

/(«) s ~ 52 P J lo « P i- < 4 ' 10 ) 

3 = 1 

Note that if each cell has the same visitation probability, so that Pj = l/N(e) for all 
j , then 7(e) = N(e) and Di reduces to Dc- In general, though, 7)/ < Dc- 

To generalize, it is convenient to define an infinite set of dimensions. Let X(t) = 
(ari(t), x 2 (t), ■ ■ ■ ,Xj(t)) denote a trajectory in a d- dimensional phase space and con- 
struct from the trajectory a sequence of points ^(0), ^f(r), X(2t), ■ ■ ■ , X(Mt) in 
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phase space, where M is large. Partition the phase space into N(e) cells of side t 
and let A/, be the number of points appearing in cell i. Define P, = MJM as the 
probability of finding a point in cell i. Then the dimensions D n are defined by 

m \ 

E^ n ]/loge. (4.11 

Note that 

A) — — lim log N(c)/ log c = Dc, (4-12) 

£—►0 v ' 

and 

/*(«) \ 

D i = Jim I Y, a log a I / log e = D v (4.13) 

That is, D q and Di reduce to the fractal and information dimensions defined earlier. 
For the total set of generalized dimensions it can be shown that D n < D m for n > m, 
where the equality holds for a completely homogeneous probability distribution, i.e. 
Pi — 1 / N (e) . Hence, the difference between dimensions of different order measures the 
degree of inhomogeneity of the attractor due to the variation of visitation frequency 
of each cell. 



The dimension D 2 is related to the correlation integral function 


1 


M 


tj=l 

(«¥i) 


(4.14) 


where is the separation of two of the points such that z i3 = [J2k=i( x k(tj) ~ 
( A ) ) 2 ] 1 ^ 2 anc l @( x ) I s the step function such that 


9(x) 



if x < 0 
otherwise. 


(4.15) 


To see this, note that P? is the probability that any two points lie in the same cell. 
This is approximately the probability that two points on the attractor are within a 
distance t of each other, which is given by eq. 4.14. Thus, the correlation dimension 
is 


D 2 (d ) - lim 

e^O 


logA M) 

logo 


(4.16) 
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4.2.1 Kolmogorov entropy 


Another useful quantity to characterize a chaotic attractor is the Kolmogorov 
entropy, which involves a concept with a strong analogy to the entropy concept in 
information theory. In information theory, we know that the maximum information 
content for a system with two possible states is 1 binary bit. It is 2 bits for a system 
with four possible states and so on. In general, the maximum information content is 
log 2 N bits for N possible states. 

Suppose a system has N states, and the ith state is known to occur with probability 
Pi. The entropy of the system is defined as 

S = -jr Pl \og 2 Pi. (4.17) 

i= 1 

Here, entropy is a measure of the amount of information necessary to determine the 
state of the system. Consider now a trajectory x(ti ) = (ac 1 (t t -), ac 2 (^t)» • • • i x n(U)) and 
partition phase space into n hypercubes of side e. Let P(ii,i 2 , • ■ • An) be the joint 
probability that the point x(ir) lies in the i;th cell. Then, from the above discussion, 


K n = - ^(*l»*2,---»*n)log 2 P(*i,*2,...,*n) 

*1 «•••.*?» 


(4.18) 


is a measure of the amount of information necessary to specify the trajectory to within 
a precision e. It follows that K n+ i — K n is the additional amount of information 
required to specify which cell x(nr + r) will fall in. The Kolmogorov entropy may be 
defined as [63] 


K 


limlim Jim^ ^(K n +i - I< n ) 

-limlim lim Y\ P(h, * 2 , • • • , *n) 

T—0 e-*0 N-> oo \NtJ ■ 

tl ,«2 *n 

^ ^§2 


(4.19) 


(4.20) 


We see that K is the average rate of the information loss. Hence, K = 0 in an ordered 
system; K = oo in a random system; and K is positive constant 0 in a chaotic 
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system. To generalize, we define the order-tf Renyi entropy as 

' 1 


K q = — lim lim lim 

r— +0 c— fO TV— foo 


Nr{q 


jT ) lo §2 E *2,- (4.21) 

'' t‘i ,12,..., i n 


with P 9 = Pexp(q - l)log 2 P. Of all the order-? quantities K q , I< 2 is the easiest 
quantity to calculate. It can be shown to be a lower bound to the Kolmogorov 
entropy [63]. 


4.2.2 Dimensions and phase space portraits from an experimental time 
series 

In the discussion of the preceding section, we made no mention of how one can 
construct a phase space portrait of a dynamic system with only time series data for one 
variable. For example, in our experiments with the PCR, the data is the temporal 
fluctuation of an intensity. According to the Takens embedding theorem [64], the 
topology of an attractor can be reconstructed in an artificial phase space derived 
from a single experimental time series. The procedure in such cases is to construct 
a pseudo phase space, for example by plotting the value of the measured variable at 
time t on one axis and successively time-delayed values on other axes. For example, 
we denote the measured quantity at time t by x(f). Then, for a two-dimensional 
phase space, we plot x(t + r) on the second axis, where r is a small interval of 
time. Construction of a d embedding dimensional phase space follows in the same 
manner by using X = (x(t),x(t + r), • • • ,x(t + dr)). The choice of r is a debated 
subject [56]. There is a certain amount of trial-and-error involved in finding the best 
r for a particular problem. If r is made too large, one loses the correlation between 
the data points, with the result that the plot is just a randomly moving wiggly line. 
One the other hand, if too small a value is chosen, the plot shows strong dependence 
among points and no structure can be resolved. 

To calculate the correlation dimension and the Kolmogrov entropy from an exper- 
imental time series data, we simply apply the calculation methods discussed in the 
preceding sections to the embedding phase spaces made up of those vector points X . 
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The second order Renyi entropy is given by 

K 2 = r -1 [logC 2 (e,<i) - log C 2 (e, d + 1)]. (4.22) 

For small e and large embedding dimension d it is found that 

C 2 (e, d ) ~ t D 2( ' d ' 1 exp(— dr/lf 2 ), (4-23) 

so that in practice K 2 can be obtained by calculating C 2 {e,d) for different values of d. 
As I< 2 is a lower bound to the Kolmogorov entropy, the condition I< 2 > 0 is sufficient 
for the existence of chaos. 

4.3 Background on spatiotemporal chaos 

4.3.1 Extended systems. Weak turbulence 

Chaos, as discussed in the preceding sections, is a temporal phenomenon as long 
as there are no spatial counterpart involved. This limitation has to be reconsidered 
if the degrees of freedom in a dynamical system have a spatial significance. Indeed, 
in nonlinear nonequilibrium states, as in fluid convection, there appear structures, 
which have been called dissipative structures by I. Prigogine [65]. Examples of these 
states having structures include: 

(1) patterns with steady structures which vary purely in space (e.g. convection of 
fluids, chemical reactions, etc.); 

(2) states with homogeneous structures which vary periodically purely in time (e.g. 
chemical oscillations, heartbeat etc.); 

(3) quasiperiodic states (similar to (2)), which vary quasiperiodically in time; 

(4) chaotic states which vary chaotically in time; 

(5) turbulent states which vary chaotically both in time and in space. 

Some physical systems may only exhibit one of these examples while others can 
have all these, such as Rayleigh- Benard convections and Taylor-Couette flows [66]. 
Regarding the nonequilibrium states involving spatial structures, one of the tasks of 
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spatiotemporal chaos is to find out how the local/global instability mechanisms are 
manifested in the physical space rather than in the phase space portraits. 

The notion of temporal chaos discussed in this thesis applies to a situation where 
the confinement effects are so strong that the degrees of freedom in a dynamical 
system one studies do not have any spatial significance. In that case, one deals only 
with the temporal aspect of the dynamics. When systems are weakly confined, i.e. the 
boundary effects are less stringent, the instability mechanisms can be then manifested 
in the structural modulations. The word pattern will then be a key to the associated 
disorder, which has now a space-time meaning. 

A much studied example of a confined system is that of a fluid layer contained 
in a small box heated from the bottom [67]. At the beginning of the experiment, 
the liquid is static. I he heat flows across the layer cell via conduction. As the 
temperature gradient reaches a critical value, a convection flow sets m. The hot and 
less dense liquid rises in the middle while the cool liquid flows down at the sides. So, 
there appear two convective rolls. As the temperature gradient is increased further, 
the rolls suddenly become unstable in a specific way. A wave (mode) starts running 
along the roll (see figure 4.5) periodically. As the temperature gradient is increased 
further, an additional new wave (mode) is observed. The temperature monitored 
locally fluctuates periodically and quasiperiodically until it becomes chaotic. This 
early work has shown how the confinement effects affect the dynamics. From the 
physical dimensions of the container, the aspect ratio T can be defined to be the 
ratio of the laterial extension / to the height which is comparable to the radii of the 
rolls, A r , i.e. T = lj\ r . When T is small, the spatial structure is stablized by the 
lateral confinement effects. In constrast, when T 1, the system is called weakly 
confined or extended, where the size of ordered states is down to a local scale only. 
The overall structure is then to better understood in terms of pattern or textures. As 
the complexity level increases, structural defects may often appear. The role played 
by such defects in the turbulent system is a currently active area of research [68]. 

If now the layer of fluid considered above is contained in a box of much larger 
aspect ratios, i.e. the lateral extension is much larger than that of the radii of the 
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Figure 4.5: Schematic diagram of the Rayleigh-Benard convection. A thin layer 
of fluid in a small box is heated from the bottom. Two convective rolls rotate as 
indicated by the arrows when the temperature gradient reaches a critical value. As 
the temperature gradient further increases, waves are excited and run along the lateral 
direction of the rolls. 
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rolls, a continuous spectrum of excited waves (modes) can be present. When many 
of them with comparable wavelength can be present independently with nearly equal 
amplitude, their mutual interferences can induce large scale modulations. These 
kinds of modulation have been associated closely with the so called the Benjamin- 
Feir instability [69]. To describe such modulations or instability mechanisms, the 
basic idea of spatiotemporal chaos is to extend the notion of an amplitude that is 
time-dependent-only, as in temporal chaos, to account also for spatial dependence. 
In other words, the amplitude of the wave has to be described by the two variables 
of time and space. In fact, such an idea of using envelop functions for the amplitudes 
of the excited waves leads to a complex Ginzburg- Landau equation to describe the 
dynamics of the dissipative structures in hydrodynamics [68]. 

4.3.2 Topological turbulence 

The last few years have witnessed very extensive explorations of the spatial be- 
havior of laser beams which may display spatial structures similar to those seen in 
hydrodynamics [70, 71, 72, 73, 74, 75] . A series of investigations of transverse ef- 
fects in laser dynamics by Lugiato and co-workers shows that, just as temporal chaos 
results from the presence of at least a few competiting frequencies in a nonlinear 
system, instability mechanisms of the spatial structures in a laser may be governed 
by the interactions and competitions among a few chosen transverse cavity modes. 
The solutions of Maxwell-Block equations including diffraction [70] display families 
of transverse patterns containing regularly distributed dark spots called phase singu- 
larities [76, 77]. Also, recently, in their analysis of laser cavities with large Fresnel 
number, P. Coullet et al. reduced the Maxwell-Block equations to a single equation 
of the form of the complex Ginzburg- Landau equation [78]. In the two-dimensional 
space of the Ginzburg- Landau equation, they found singular solutions in the form of 
spirals. As the spiral waves bear some analogy with the vortices in hydrodynamics, 
they are now named optical vortices. In the stable case, these solutions persist in- 
definitely owing to their character of topological defect. In the unstable case, they 
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nucleate spontaneously as temporary singularities of the moving wavefront. The spa- 
tiotemporal chaotic regime that sets in has been called defect-mediated turbulence 
or topological turbulence [79, 80]. Defect- mediated turbulence is currently an active 
area of investigation across many disciplines of nonlinear science. The role of de- 
fects in mediating turbulence in nonequilibrium systems with large aspect ratios is 
an important matter of interest. Therefore, some knowledge of defects [81] becomes 
necessary. Before discussing some of the structures and characteristics of defects, let 
us close this brief introduction on spatiotemporal chaos with some mention of the 
methodologies commonly used for its analysis. 

4.3.3 Methodologies of spatiotemporal complexity 

At present, there are still no well-established methods for studying spatiotemporal 
chaos, in particular the fully developed turbulence. However, we may see emerging 
some important paradigms or methods for the study of the spatiotemporal dynamics. 
We now describe briefly some of these methods. 

One of the oldest examples of a dynamical system which may evolve both in space 
and time are cellular automata. Cellular automata were invented in the late 40s by J. 
von Neumann to mimic the behavior of complex, spatially extended structures. W. 
Burke gave a good overview of the work of von Neumann [82]. In cellular automata, 
space is divided into discrete small units called cells or sites. The cells take on 
numerical values. Consider a one-dimensional cellular automaton consisting of a row 
of cells. Each cell is assigned a set of initial numbers. A set of rules specifies how 
these numbers are to be changed at every clock step. Suppose that the automaton 
initially is filled with 0’s except for a single cell which is occupied by a 1: 

... 010000000 ... 

And suppose that the rule states that the number in each cell is to be replaced by 
the sum of itself and its left neighbor. Thus, after one clock step, the state of the 
automaton will be as follows: 

... 011000000 ... 
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Another clock step later the state will be 

... 012100000 ... 

followed by 

...013310000... 

...014641000... 


etc. In fact, this cellular automaton is one of the early von Neumann machines used 
to compute the binominal coefficient of ( a + b ) n . Cellular automata can be applied to 
one, two and higher dimensional problems. The fixed rules which specify the change 
of the cells’ values carry out local interactions between neighboring cells, reflecting 
the dynamics of the system in question. Cellular automata have been used to model 
the complexity in biological systems from the level of cell activity to the levels of 
clusters of cells and populations of organisms. In chemistry, cellular automata have 
been used to model kinetics of molecular systems and crystal growth. In physics, 
they have been used to study dynamical systems such as the clustering of galaxies. 

Another model for studying spatiotemporal chaos in spatially extended systems 
are the coupled map lattices introduced by K. Kaneko [83]. Coupled map lattices 
are dynamical systems discrete in both time and space. Although their modeling 
structure looks similar to that of cellular automata, the change of the state of each 
site is done by a local mapping function rather than by a set of rules. One example 
of coupled map lattices has the form 

*n+i(0 = (1 - D)f(x n (i)) + D[f(x n (i + 1) + f(x n (i - 1))], (4.24) 

where n represents the time step, i represents a lattice site, and D represents the 
diffusive coupling between nearest neighbors. The function /(re) is chosen to be 
ax{\— x) (i.e. coupled logistic lattice) or x + A sin(27rx) + c (i.e. coupled circle lattice). 
It has been shown that transitions to turbulence via spatiotemporal intermittency 
arise in these systems [84], At the same time, studies of spatiotemporal chaos by 
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coupled map lattices have been growing rapidly because calculations are efficient and 
they lend themselves to the generalization of some of the quantitative methods used to 
characterize temporal chaos (e.g. Lyapunov spectra and Kolmogorov entropy density 
). This approach is also one of the very few known approaches providing means for 
the diagnosis of experimental data of spatially extended systems [84], 

A third possible approach to the study of spatiotemporal dynamics is called the 
modal expansion of fields with transverse structures, as used in laser systems by 
Lugiato and co-workers. As is well known, spatiotemporal problems of continuous 
variables are mostly described by partial differential equations (PDEs). The solutions 
often rely on numerical methods. Sometimes, PDEs can be formidably demanding in 
computing power. Lugiato and co-workers demonstrated an alternative way of solving 
the PDEs of lasers by expanding the dynamical quantities in a series of normal modes. 
The PDEs are then reduced to a set of coupled ordinary differential equations (ODEs) 
in the limit of the mean-field approximation. The basic idea of this approach is to 
represent the envelop functions of the fields’ amplitudes by the space-dependent-only 
normal modes and the corresponding time-dependent-only coefficient if it is validly 
allowed. This method is highly efficient and yet the essential dynamics of the spatial 
structures can still be captured. 

A last approach to a description of the turbulent state in a nonequilibrium system 
is to characterize it by the presence of topological defects. As mentioned before, this 
approach of topological turbulence was first introduced by P. Coullet et al. for laser 
systems [80]. They used numerical solutions of the Ginzburg-Landau equation to 
illustrate a mechanism of defect formation, annihilation and dynamics in dissipative 
structures. Another possible theoretical approach to describe defect formation and 
turbulence in large aspect ratio systems has recently been proposed by I. Procac- 
cia [85]. In the study of nematic electroconvection, the ordered state of the convec- 
tive rolls before the appearance of defects is well described by the envelope equations. 
However, how the modulations are affected by the defects is not clear. Procaccia 
concluded in his study of nematics that a field theory might be necessary to describe 
the dynamics and the interaction forces between defects. On the experimental side, 
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several examples have already demonstrated defect formation in nonequilibrium sys- 
tems of large aspect ratios, as in surface waves, in thermal convection, and in the 
transverse patterns of lasers and nonlinear optical systems [56, 86]. In fact, the dy- 
namics that we observe in the PCR above threshold also provide evidence for such 
defect-mediated transitions [87, 88]. This will be demonstrated in the following chap- 
ters. But we first introduce in the next section some of the concepts related to the 
definitions of defects and in particular the phase defects or wavefront dislocations 
observed in electromagnetic waves. 


4.4 Wavefront dislocations/phase defects 

In this concluding section, we will define the singularities found in electromagnetic 
waves and the corresponding terminology used later on in our discussions of the PCR’s 
dynamics. The discussion given here follows closely that of a paper by M. Berry [89]. 

Wavefront dislocations or phase defects were first discovered by Nye and Berry 
in their study of ultrasound reflected from a rough surface [90]. Wave dislocations 
appear ubiquitously in optical fields, acoustics, and quantum mechanics as well as 
water waves [81, 91]. Wavefront dislocations are just one kind of singulary that can 
appear in scalar waves (for vector or tensor waves, see [92]). Ray caustics are another 
kind of wave singularity with a different length scale. On a caustic, rays of light (in a 
geometrical optics sense) are all focused to a point and thus the intensity is infinite. 
On a phase defect, the intensity is zero. To be observable, the phase defects of a 
wavefront must be explored on the wavelength’s scale, for example by interferometry. 
Caustic singularities, in constrast, are visible to the naked eye (e.g. the reflection of 
sunlight from a wavy water surface). 

As is well known, singularities are present in many complicated systems, such as 
solids, liquid crystals, and hydrodynamics as well as in nonlinear optics. When con- 
sidered as geometric structures, singularities are forms rather than materials. Waves 
are also forms. Therefore, dislocations in wavefronts are forms of forms. Despite their 
“double” abstract nature, they are directly observable phenomena. 
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Another important nature of wavefront dislocations is their singular character. 
Consider a travelling wave ip(r,t) with complex amplitude p(r, t) and phase x( r >0 
defined by 

ip = pexp(ix)- (4.25) 


For the waves in which we are interested, ip is a smooth single-valued function of 
position and time. When traversing a loop in the wavefront, \ ma y change by 2m7r, 
where m is an integer. Let us shrink the path to a very small loop while the amount 
of variation in x does not change. Thus, a singularity is enclosed because x is varying 
infinitely fast. And yet the smoothness of ip has to be preserved. It implies that this 
can happen only where ip = 0, i.e. x is indeterminate. As ip is null, both R eip and 
Im ip vanish (i.e. R eip = Imip = 0). Thus, phase singularities are points in planes and 
lines in space, analogous to crystal dislocations in solids. 

Before we illustrate the possible structural forms of the wavefront dislocations by 
examples, we would like to ask how they happen. Since the singularities happen only 
where ip = 0, i.e. the intensity p 2 vanishes, they might be confused with the dark 
fringes often discussed in elementary treatments of interference. The dark fringes 
are, however, conceived as surfaces rather than lines. Wavefront singularities result 
from the destructive interferences of the different rays contributing to ip at a crossing 
region of waves scattered from different objects. The crossing can be a point in a 
plane or a line in space. 


Now we will illustrate some of the possible structural forms of the dislocations of 
a wave ip satisfying the wave equation in a uniform static isotropic medium: 

1 d 2 ip 


V 2 ip = 


c 2 dt 2 


(4.26) 


Suppose that ip travels only in the z-direction, i.e. ip = pexp(z£), £ = z — ct. The first 
case considered is a dislocation moving rigidly with the host wavefront. The wave for 
such a case must satisfy 

d 2 ip d 2 ip 


ip = ip(x,y, 0, 


+ 


- 0 . 


(4.27) 


dx 2 dy 2 

Consider as an example, ip = A(kx + i0k£)e ik t, where A and j3 are real constants. 
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Figure 4.6 shows the contours of (, i.e. wavefronts of ip moving upwards. The 
dislocation is located at E. In crystallographic terms, this is a pure edge dislocation. 
The next example of dislocation is the pure screw dislocation, which for strength s 
has a wave function 

ip = A(kx ± iky) s e lk (4.28) 

Figure 4.7 shows the typical plot of pure screw dislocation with strength s = 2. 
The sign ± determines the helicity of the screw and the charge s its pitch. Negative 
charge of a screw dislocation means that the phase increases counterclockwise around 
the dislocation. There is also a mixed type of dislocation, i.e. edge-screw dislocation 
with waves of the form 

ip = A(kx + i(3(k£ cos 8 — ky sin 6))e ,fc ^, (4.29) 

for example, (5 = 0 and 6 = n/2 correspond to pure edge and pure screw respectively. 

Now consider the modulation of the wave in equation 4.25 with quadratic terms in 
x, y and which also satisfy equation 4.27. In this way, it would have curved dislocation 
lines that move rigidly with the host wave. An example is ip = [(x — iy) 2 + i(3k(,]e' k t 
whose dislocation lines have the form of two parabolas orthogonal to each other, as 
shown in figure 4.8a. 

We now consider examples showing dislocations moving relative to the wavefront 
as the host wave propagates. A simple example of these is 

tP = A[akx + k 2 x 2 + i(0kt + kz)}e' k t, (4.30) 

where a is real constant. This has two edge dislocations at x = 0 and x = cxjk^ 
moving parallel to the propagating axis (see figure 4.8b). In a crystallographic term, 
such motion is glide. Another kind of movement is made possible if the term in £ is 
made quadratic, i.e. if 

ip = A[akx + k 2 x 2 + i((3k 2 t; 2 + kz)]e' k (4-31) 

Now two pairs of edge dislocations are located at x = 0, and x = — a/k . If /? > 0, 
members of each pair approach, collide and annihilate while if /? < 0 they appear and 
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Figure 4.7: Wavefronts of a travelling wave containing a pure screw dislocation. This 
helical wavefront with the dislocation along the axis of propagation rotates in clock- 
wise. The direction of rotation determines the sign or charge of the screw dislocation. 
Clockwise rotation means a positive charge and vice versa. 
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Figure 4.8: Examples of dislocations’ motion in a travelling wave. The wave travels 
upwards (see text for details). 
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Figure 4.9: More examples of dislocations’ movements in a travelling wave. The wave 
travels upwards (see text for details). 
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then glide apart (see figure 4.9a). If a slight change of form is made in equation 4.31, 
for example, if it becomes 

i/> = A[k 2 x 2 + i((3 — i^k 2 ^ 2 + ikz\e' k ^, (4.32) 

where 7 is a real constant, two pairs of edge dislocations approach each other along 
parabolic trajectories and annihilate when they meet (see figure 4.9b). Such motion, 
called climb, corresponds to the spontaneous recovery of a tear in a wavefront. Above 
are a few of many possible structural forms and movements that the wavefront dis- 
location or phase defects can have. Many other interesting examples can be found in 
references [90, 93, 89]. 

In conclusion, we have presented a brief introduction of nonlinear dynamics and 
spatiotemporal as well as temporal chaos. Of particular interest are the definitions 
of topological turbulence and the role played by the wavefront dislocations or phase 
defects in mediating the dynamics. Some aspects of these phenomena will be illus- 
trated in the study of the PCR’s dynamics above threshold, described in the following 
chapters. 
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Chapter 5 


Experiments: spatiotemporal 
dynamics in PCRs 


In this and the following chapter we present the results of an experimental study 
of spatiotemporal structures in the PCR’s beams above threshold. In this chapter, 
the emergence of spatiotemporal structures as the degree of transverse confinement is 
varied by changing the Fresnel number will be described in detail. A theoretical model 
of the spatiotemporal dynamics of the PCR will be presented in the next chapter. 

This chapter is organized as follows: In section 5.1, the experimental arrangement 
for the study of the spatiotemporal dynamics of the PCR will be described. The 
experimental results and their analyses are presented in the following five sections. 
The dynamics of the local intensity measured at one point in the beam will be analysed 
in section 5.2 with the aid of experimental time series, power spectra and phase 
space portraits. An attempt at quantifying the chaotic dynamics by calculating the 
correlation dimension and the Kolmogorov entropy will be presented in section 5.3. 
The spatial instabilities of structures across the PCR’s beams will be discussed in 
sections 5.4 and 5.5, where it is indicated that the spatiotemporal dynamics observed 
has some the characteristics of a topological or defect-mediated turbulence. Finally, 
in section 5.6, it is shown that the superposition of a few simple empty cavity TEM 
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modes with different frequencies have a dynamical behavior compatible with the 
observed phenomena. Additional discussion and conclusions are presented in the 
concluding section 5.7. 


5.1 Experimental setup 


The linear phase-conjugate resonator under study is sketched schematically in 
figure 5.1. The resonator consists of a flat dielectric mirror M (95 % reflectance) and 
a phase-conjugate mirror. The phase- conjugate mirror is a single crystal of BaTi0 3 
5mm x 5mm x 7mm in size. It is pumped externally by two counterpropagating beams 
from a single-mode Argon ion laser tuned to 514nm. The beam diameter is about 
1.5mm. The intensity and polarization direction of each pump beam are controlled 
independently by assemblies consisting of a half-wave plate and a crystal polarizer. 
In all experiments reported in this chapter, both beams’ polarization directions were 
extraordinary, i.e. in the same plane as the c-axis of the crystal, which lies in the 
plane of the figure. 


The resonator contains two apertures, Hi and H 2 . The aperture H 2 is located as 
close as possible to the mirror M and in the front focal plane of a 16 cm focal length 
lens. The aperture Hi is in the back focal plane of that lens, about 3 cm away from 
the phase-conjugate mirror. The degree of transverse confinement of the system can 
be varied by changing the apertures sizes. A quantitaive measure of this confinement 
is given by the resonator Fresnel number defined as 


F = 


D\D 2 

"Tp 


(5.1) 


where D x and D 2 are the diameters of the apertures Hi and H 2 respectively, A is the 
wavelength and / is the focal length of the lens. Di is small enough to ensure that 
the interaction region in the crystal is smaller than the area of pump overlap. The 
Fresnel number has a concrete physical meaning [94]. A simple argument based on 
scalar diffraction theory shows that the transverse dimension of a speckle near the 
aperture Hi is given by A f / D 2 . F is, therefore, a measure of the number of speckles 


98 



CHAPTER 5. EXPERIMENTS: SPATIOTEMPORAL DYNAMICS IN PCRS 


D3 



Figure 5.1: Schematic diagram of the phase-conjugate resonator setup: A/2-P are 
half-wave plate - polarizer assemblies used to vary the intensity of the beams indi- 
vidually. All beams’ polarizations are extraordinary. The two apertures H 1; H 2 are 
one focal length away from the lens (f=16cm). The Ar + laser operates in a single 
mode at 514nm. The phase-conjugate mirror is a single crystal of BaTi0 3 . Dl, D2 
are photomultiplier tubes measuring the local intensity at two different locations; D3 
monitors the laser output and the CCD camera is used to capture the dynamics of 
the intensity distribution. 
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along any transverse diameter in the resonator field. In all the experiments reported 
here, except for the measurements made with F < 3 , the aperture Hi was fixed with 
D\ = 0.69mm. 

The Fresnel number is the control parameter used in the experiments. The Bragg 
detuning achieved by misaligning one of the two pumps is another significant param- 
eter, which will be considered in the next chapter. Both appear to be the parameters 
that most directly affect the spatiotemporal dynamics of the system. Two other pa- 
rameters have been observed to affect the dynamics as well, but to a much lesser 
degree. The first is the total intensity I a incident on the crystal. The response time 
associated with the photorefractive effect is roughly inversely proportional to the in- 
tensity in the range used ( lOmW/mm 2 ~ 100 mW/mm 2 ). For our crystal, the 
response time was of the order of Is mm 2 /mW. The second factor influencing the 
dynamics is the pump ratio, defined as R p — I 2 / hi where h and I\ are pump 2 and 
pump 1 intensities, respectively (see figure 5.1). The pump ratio can be used to vary 
the departure of the cavity from threshold, which affects the buildup and decay rates 
of the cavity [34]. It has been seen in chapter 2 that the decay time diverges for 
the values of the pump ratio at which the cavity reaches threshold for self-oscillation. 
From chapter 3 (fig. 3.10), it is seen that the cavity is above threshold for 0.1< R p < 
1 . The departure from threshold can then be varied by choosing a pump ratio within 
that range. With a better alignment of the crystal, and thus a higher gain coefficient, 
that range could be extended to R p < 2 for the experiments described in this chapter. 

With a tight confinement (F = 1.8), giving a stationary output, we measured the 
intensity and the buildup times of the phase conjugate beam / 3 against the pump ratio 
for various values of the total pump power. The goal was to determine what values of 
these two parameters gave rates of change of the observed phenomena compatible with 
visual observation and standard video recording. The results are shown in figure 5.2. 
For all the experiments reported here, the pump ratio was fixed at 1 and the total 
power of the two pumps was kept at 60mW. Other parameters such as alignments, 
four- wave-mixing geometry, etc., were fixed for all the measurements. 

The field distribution in aperture Hi was imaged (with lenses not shown in Fig. 5.1) 
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Figure 5.2: (a) Phase conjugate intensity / 3 versus the pump ratio for various values 
of the total pump power, (b) Buildup time of / 3 versus pump ratio for various values 
of the total pump power 
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through the crystal onto three output planes. The local intensities at two different 
locations in Hj were measured through ports 1 and 2 by two photomultipliers D1 and 
D2. The points at which the intensities are measured were selected by two pinholes. 
Both photomultiplier currents were sent to current-to-voltage amplifiers and to analog 
filters before being digitized and recorded by a PC based system. 

The transverse intensity distribution and its dynamics was recorded via port 3, 
either by a motorized 35mm camera or a CCD camera connected to a video-tape 
recorder. The same port and the same modes of recording were also used to record 
interferograms of the wavefront. To do this, the field in the image plane at port 3 
was made to interfere with a tilted plane wave from the same laser. 

The cavity output through the phase-conjugate mirror was chosen because of its 
higher intensity. It is amplified in the crystal. The transverse beam profile and its 
dynamics were found to be the same in the entire volume of the cavity. This was 
verified by extracting the left and right propagating cavity beams with the aid of a 
low loss pellicle beam splitter, and comparing the profiles of these beams with that of 
the output used for the measurments. The phase conjugate beam exiting the cavity 
via the conventional mirror had also a similar transverse profile and exhibited the 
same dynamical behavior. 

An additional detector D3 was used to monitor the laser output and check for any 
spurious correlation between the measured time series and the laser output fluctua- 
tions. The output was found to be stable with occasional but very rare mode hops. 
Since such events would affect the cavity dynamics unpredictably, the data would be 
trashed when a mode hop occurred during the recording of a time series. 

5.2 Time series, Power spectra and Phase space portraits 

As we have seen in section 4.1, the state of a dynamical system changes with the 
value of the control parameters, and, in some cases, the system may become chaotic. 
Time series, power spectra and phase space portraits are useful tools for examining 
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the dynamics of a nonlinear system. 

In our experiments, time series of the local intensity containing 8,000 data points 
each were recorded for various values of the Fresnel number. Some representative 
examples are shown in figure 5.3. A number of characteristic transitions from order 
to chaos can be identified in this figure. 

For tight confinements^ ^ 2), the irradiance distribution in aperture Hj observed 
through the phase-conjugate mirror is uniform and the local intensity is stationary. 
A first bifurcation occurs near 2 < F < 2.2, leading to stable periodic oscillations. As 
the confinement is relaxed further (2.2 < F < 3.5), more bifurcations occur, leading 
to more complex but still periodic motions. An example is shown for F = 3.2. 
Eventually, these higher order bifurcations lead the system to a chaotic state near 
F ~ 4. At higher Fresnel numbers (F > 4), the nature of the dynamics appears to 
change. At F = 5.8 for example, a quasiperiodic state interrupted by chaotic bursts 
is being observed. Finally, at F > 7, a new chaotic state sets in. 

Time series of 8,000 points may be considered as short, especially near bifurca- 
tion or threshold points where critical slowing down might take place. Longer time 
series, however, would have had a much higher probability of corruption from laser 
instability, instrumental drift, and mechanical relaxation. 

The power spectra of the time series of figure 5.3 are shown in figure 5.4. These 
spectra reveal information complementing that contained in the time series. 

At F — 2.2, the motion is clearly periodic with a fundamental frequency f!=0.106Hz 
and its higher harmonics. At F — 3.2, the spectrum contains, in addition to the fun- 
damental frequency f!=0.132Hz, all the multiples of the subharmonic fj/3. This is 
the spectrum of a period three motion. The corresponding time series indeed shows 
three peaks per period. 

As F is increased further and the transverse confinement is relaxed, additional 
subharmonics appear. At F ~ 4 however, a dramatic change occurs. The main 
feature of the spectrum is now a broadband continuum, a clear sign of an irregular 
signal. At low frequencies, near 0.055Hz where some new feature seems to have 
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Figure 5.3: Output intensity versus time at a point within aperture H x , for different 
values of the Fresnel number F. The total power of the two pumps is 60mW and the 
pump ratio is 1. The time series shows (a) a stationary output, (b) stable oscillations, 
(c) the emergence of subharmonics, (d) the development of chaotic oscillations, (e) 
unstable oscillatory motions with intermittant bursts, and (f) another chaotic motion. 
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Figure 5.4: Normalized power spectra calculated from the data of Fig 5.3 showing 
(a) a fundamental frequency and its higher harmonics, (b) a fundamental frequency 
and its subharmonics multiples of 1/3, (c) a broadband spectrum with the emergence 
of two frequency peaks, and (d) two broadened frequency peaks and their linear 
combinations. 
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appeared in the spectrum, the background is nearly three orders of magnitude above 
the noise floor. One can thus rule out random noise external to the system and safely 
attribute the appearance of the continuum to the intrinsic dynamics of the system. 
At higher frequencies, the power is seen to drop linearly in the semilog plot. This 
indicates an exponential decay of the power spectrum. Such a spectral feature has 
been recognised as a signature of deterministic chaos [95]. 

At F = 5.8, the spectral energy is mostly confined near two different frequencies, 
f!=0.128Hz and f 2 =0.055Hz, and their linear combinations. The continuous back- 
ground has not entirely disappeared, however. It is interesting to note that these two 
frequencies can already be seen to emerge from the background at F = 4.1. This 
identification of two frequencies, one of them being the frequency of the original limit 
cycles and another, apparently unrelated, suggests a route to chaos following the 
Ruelle-Takens-Newhouse scenario. 

As seen in section 4.1.2, the pseudo phase space portraits provide additional infor- 
mation on the dynamical state of the system. According to the Takens’s embedding 
theorem, a pseudo phase space portrait having the same geometrical properties as 
the original one can be reconstructed from a single time series of experimental data 
(section 4.2.2). Some examples of pseudo phase space portraits are shown in figure 5.5. 

At F = 1.8, the phase space portrait is a simple fixed point. A limit cycle is 
the portrait at F = 2.2 as shown in figure 5.5a. The width of the attractor may 
be attributed to stochastic noise. At F = 3.2, a stable period three cycle appears 
(figure 5.5b). The phase space portrait at F = 4.1 is an irregular distribution of 
points filling up the phase space more or less uniformly (figure 5.5c). The phase space 
portrait at F = 5.8 (figure 5.5d) exhibits some very interesting features indicative of 
what appears to be an unstable limit cycle surrounded by a cloud of scattered points. 
As seen in figure 5.6, this phase space portrait corresponds to a time series which 
is nearly periodic with occasional interruptions by irregular bursts. This behavior is 
compatible with the Ruelle-Taken-Newhouse scenario in which additional frequencies, 
unrelated to that of the original limit cycle appear in the spectrum. In this scenario, 
the bursts are interpreted as a manifestation of type III intermittency [60]. Indeed 
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Figure 5.5: Phase space portraits of time series containing N=8000 data points. I(i) 
is plotted versus I{i + n), with a time delay given by At = nr. r is the sampling 
interval. One can identify (a)a limit cycle, (b)a stable period 3 motion, (c)an irregular 
filling of phase space; and (d)an unstable limit cycle smeared by a cloud of irregularly 
distributed points. 
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some of the intermittent bursts appear to show an increase and then a sudden drop 
of the amplitude of a subharmonic, which is a signature of type III intermittency. 

5.3 Correlation dimension and Kolmogorov entropy 

We have seen in section 4.1 that a number of quantitative methods have been es- 
tablished for studing dynamical systems and distinguishing deterministic chaos from 
random noise. In the following, we analyse quantitatively the nature of the chaotic 
motion observed with a Fresnel number F = 4.1 by calculating the correlation dimen- 
sion and the corresponding entropy, using the corresponding time series, a section of 
which is shown in figure 5.3d. 

The Grassberger-Procaccia algorithm was used to calculate the correlation di- 
mension of the experimental chaotic time series [96]. This algorithm is based on the 
Takens embedding theorem according to which N vectors X(t n ) = [X(t n ),X(t n + 
At), . . . ,X(t n + (d — l)Af)] of length d ( the embedding dimension) can be extracted 
from the time series and used to reconstruct the attractor in a pseudo d-dimensional 
phase space. The delay At = ir, where r is the sampling interval, is chosen such that 
the choice of r is made according to the remarks explained in section 4.2.2. 

The correlation function or second order correlation integral is then calculated as 
(see section 4.2): 

C2 M) = 1^ £ 0(e-*m.»), (5.2) 

m,n= 1 
d - 1 

where 6 is a step function and z m<n — y^|A n+ , — X m+i \ 2 is the Euclidean distance 

: = 0 

between the pair of vectors m, n and d is the embedding dimension. The correlation 
integral thus measures the number of pairs of vectors with a distance smaller than e. 
From equation 4.23, the correlation integral scales with the embedding dimension as 

C 2 (t,d) cx e I/ exp(— drK 2 ). (5.3) 

For large enough d and small e, the correlation exponent v becomes the correlation 
dimension D 2 . K 2 is the second order Renyi entropy. These two quantities are lower 
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bounds to the fractal (Hausdorff ) dimension D 0 of the attractor and to the information 
dimension D x (D 0 > D\ > D 2 ), and to the Kolmogorov entropy K\ {K x > I< 2 ), 
respectively (see section 4.1). 

Figure 5.7 shows log-log plots of the correlation integral for values of the embedding 
dimension ranging from 7 to 27. The slope of the curves, which is the correlation 
exponent according to eq. 5.2, appears to converge toward what is believed to be the 
value of the correlation dimension D 2 . The local slope, calculated using a seven point 
local average of the slope of the data of figure 5.7, is shown in figure 5.8 for different 
values of d. Although the plateau, which should indicate the value of D 2 , is not very 
broad, it does seem to converge, for d > 20, toward a value approximately equal to 
5.2, as shown in figure 5.9. 

According to eq. 5.3, a lower bound to the Kolmogorov entropy can also be ex- 
tracted from the data of figure 5.7 by calculating the rate at which the correlation 
integral drops with increasing values of the embedding dimension. In section 4.2.1 
we discussed briefly that the Kolmogorov entropy measures the rate of information 
flow in the dynamical system and is an indication of the system’s sensitivity to initial 
conditions. In an ordered system, it is expected to be zero, while in a stochastic 
system, it should be infinite. For a deterministic chaos, its value would be finite and 
positive. Figure 5.10 shows a plot of K 2 versus d. Again, for d > 20, the entropy 
appears to converge toward a value approximately equal to 0.16 s -1 . 

Although the relatively large value found for the correlation dimension may pos- 
sibly put in doubt the validity of the algorithm used to calculate it, both the finite 
positive lower bound to the Kolmogrov entropy and the non-integer value of the 
correlation dimension seem to provide some evidence for the existence of a strange 
attractor for the time series at F = 4.1. 
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log(e) 


Figure 5.8: Local slope i/(e) of the correlation integral of Fig 5.7. versus loge. The 
slopes are averaged over 7 local points. The plateau indicates the range of linear 
scaling for each d. The slope in the region of the plateau converges toward the 
correlation dimension D 2 - 
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Figure 5.9: Plot of the correlation dimension for each embedding dimension d of 
Fig 5.7. The correlation dimension appears to converge to a value D 2 ~ 5.2. 
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Figure 5.10: Value of the entropy K 2 ,d calculated from the data of Fig 5.7, versus the 
embedding dimension d, averaged over the linear scaling region in e. The asymptotic 
value I< 2 ~ 0.16s -1 , is a lower bound to the Kolmogorov entropy. 
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5.4 Wavefront dislocations 

The transverse irradiance distribution within the aperture Hi as seen through 
the PCM was recorded as a function of time. As the transverse confinement of the 
system is relaxed, the distribution of structures across the PCRs beams becomes 
more complex, containing more and more bright lobes and dark lines around them. 
At the same time, the motion of the structures also increases in complexity, changing 
from a simple periodic motion to an apparently chaotic dance. 

Interferograms of the wavefront were recorded to reveal the presence of wavefront 
dislocations or vortices [87] (also see section 4.4). In the interferograms, the defects 
should appear as a discontinuity in the interference fringes and thus should be easily 
identifiable. The number of defects and the complexity of their motions in time were 
found to be strongly correlated to the temporal dynamics of the time series described 
in section 5.2 as well as to the spatial decorrelation of the pattern. This suggests that 
the topological defects may indeed act as the catalyst for the observed spatiotemporal 
dynamics and induce the observed turbulence in a manner similar to that of the so- 
called topological turbulence (see section 4.3.2). In this section, we will describe how 
the transverse pattern changes when the Fresnel number, and thus the transverse 
confinement of the cavity, is varied. 

After a first bifurcation at 2 < F < 2.2, the transverse pattern varies in time 
in a simple and periodic way. Some snapshots of the irradiance distribution, taken 
at a half second interval between frames, are shown in figure 5.11. The irradiance 
distribution is at first uniform, both in amplitude and phase. This lasts for about 
eight seconds, after which time a folding of the phase appears at the center of the 
pattern. This is manifested by a darkening of the area. Very rapidly, a pair of defects 
of opposite charge develops from this fold at the center of the pattern. The two 
defects split and repel each other, swiftly moving away along diametrically opposite 
directions. They finally disappear at the boundary. The presence of the defects lasts 
only about two seconds or so. The wavefront then quickly recovers and remains stable 
for a period of time. Then, the sequence repeats itself. Figure 5.12 is an enlarged 
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Figure 5.11: Sequence of snapshots of the irradiance distribution in aperture Hi with 
F = 2.2. The time interval between sanpshots is 0.5 second. 
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Figure 5.12: Instantaneous interferogram with F = 2.2 showing a pair of defects of 
opposite charge. 
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snapshot of an interferogram obtained by interfering the output distribution with a 
tilted plane wave and showing the pair of defects of opposite charge moving away 
from each other and along a horizontal diameter. 

As the transverse confinement of the system is relaxed, more complex behaviors 
are observed. At F = 3.6, for example, the motion of the transverse pattern is more 
complex but still periodic. This explains the appearance of the subharmonics in the 
time series of the local intensity (Fig. 5.3). A sequence of snapshots taken at half 
second intervals is shown in figure 5.13. Figure 5.14 is an enlarged interferogram 
showing three pairs of defects of opposite charge. 

A number of characteristic motions of the defects were observed for 2.2 < F < 5. 
Some of them are illustrated schematically in figure 5.15. Figure 5.15a corresponds 
to the motion described in figure 5.11. In fig. 5.15b, the two defects of a pair are seen 
to spiral outward, away from each other, and then inward toward each other, along 
trajectories of alternating handedness. Figure 5.15c describes a motion similar to that 
of figure 5.13. Pairs of defects are generated at two diametrically opposite locations. 
The two inner defects spiral toward each other while the other two spiral outward 
and die at the boundary. The two inner defects may rebound, repel each other and 
accelerate outward along a diameter. The sequence then repeats itself with other 
pairs of defects originating at conjugate locations and with trajectories of opposite 
handedness. Figure 5.15d shows another dance of similar nature, where several pairs 
of defects are generated, move, rotate, collide, rebound or annihilate each other. 

The motions described in the preceding paragraph may appear in sequence, when 
the dynamics is periodic. They may also appear in an apparently random order. This 
is the case for F = 4.1, which, as seen in section 5.2, corresponds to a chaotic time 
series. Figure 5.16 shows some snapshots (0.5 second interval) of the nonperiodic 
motion at F = 4.1. Figure 5.17 is an enlarged interferogram showing the presence of 
four pairs of defects. 

All the patterns with low enough Fresnel numbers have a center of symmetry. We 
believe that this is a result of the particular geometry of the cavity and of the nature 


118 



CHAPTER 5. EXPERIMENTS: SPATIOTEMPORAL DYNAMICS IN PCRS 



Figure 5.13: Sequence of snapshots with F — 3.6 (time interval = 0.5s) revealing a 
more complicated but still periodic motion. 
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Figure 5.14: Instantaneous interferogram with F = 3.6 showing three pairs of defects. 
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Figure 5.15: Sketches of some typical motions of the defects observed with 2 < F < 5. 
See text for details. 
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Figure 5.16: Sequence of snapshots with F = 4.1 (time interval —0.5s) revealing a 
nonperiodic motion (horizontal scale compressed). 
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Figure 5.17: Instantaneous interferogram with F = 4.1 showing four pairs of defects. 
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of the PCM. Indeed, one round trip in the cavity involves two operations: phase 
conjugation and imaging with unit magnification and inversion. Thus, if a certain 
distribution /(f) exists in the cavity, another /*(— f) should also be present, at least 
in steady state. 

In the range 5 < F < 6, drastic changes in the dynamics of the local intensity 
time series were observed. Similar changes are also found in the behavior of the 
defects. At F - 5.8, for example, the spatial distribution shown in the sequence of 
figure 5.18 is rather complex and indeed shows more structure than the pattern at 
p — 4.1 (Fig. 5.17), which corresponds to a chaotic motion. However, not all the dark 
lines in figure 5.18 contains true zeros. Thus, when F = 5.8, the number of defects 
is actually smaller on average and the spatiotemporal complexity is lower than when 
F = 4.1. The motion for interferograms at F - 5.8 is also simpler in appearance. 
Figure 5.19 shows an instantaneous interferogram at F = 5.8, where only one pair 
of defects is visible. When the transverse confinement is further relaxed, the motion 
becomes chaotic again (Fig. 5.3f). Figure 5.20 shows an interferogram obtained with 
^ Fresnel number F — — 8.6 and in which about eight pairs of defects are seen as they 
execute a chaotic dance. 


5.5 Defect-mediated turbulence 

It appears, from the observations described in figures 5.11 to 5.20 and the corre- 
sponding results of section 5.2, that the presence and the motion of singular points in 
the wavefront ( the phase defects or vortices ), where both the real and the imaginary 
parts of the field vanish and the phase is indeterminate (see section 4.4), may be 
responsible for the loss of spatiotemporal coherence. According to the Coullet conjec- 
ture, this in turn may induce turbulence in weakly coupled systems [79, 78, 80, 97]. 
In order to further illustrate the correlation between the number of defects and the 
spatiotemporal dynamical behavior of the system, the following two measurements 
were made. 

The temporal evolution of the interferograms obtained with increasing Fresnel 
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Figure 5.18: Sequence of snapshots with F = 5.8 (time interval = 0.5s) revealing a 
more complicated pattern but executing a simpler motion 
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Figure 5.19: Instantaneous interferogram with F = 5.8, showing only one pair of 
defects. 
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Figure 5.20: Instantaneous interferogram with F = 8.6 showing at least 8 pairs of 
defects executing a chaotic dance. 
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Figure 5.21: Maximum number N of defects observed in the transverse pattern at any 
given time versus Fresnel number F. The expected scaling law N oc F 2 is interrupted 
in the region 4 < F < 6. 
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numbers were video recorded with a standard CCD camera and examined frame by 
frame to determine the maximum number of defects appearing in the pattern at any 
instant of time. The result is shown in figure 5.21. There are two regions, 2<F <4 
and F > 6, where the maximum number of defects appears to grow as the square of 
the Fresnel number. This dependence is what is expected if, using simple diffraction 
arguments, F 2 is interpreted as being a measure of the maximum number of speckles 
( or of bright lobes ) within a cross section of the beam. The number of true zeros 
between the bright speckles is expected to be of the same order as the number of 
speckles, at least for large numbers [91, 87]. 

In addition, the temporal fluctuations of the local intensity were recorded at two 
different locations in aperture Hj, and a spatial correlation index [98] was calculated. 

The spatial correlation function is defined as 


Ctfu^T) 


N 

£['(ri,n) - /(ri)][/(r 2 ,n + f) - J(r 2 )] 


n= 1 


S(r,)S(f 2 ) 


(5.4) 


= |Z,W r '. n ) - '(r,)]'| , (5.5) 

where rq, r 2 are the two locations in the wavefront where the local intensity fluctua- 
tions are measured, I(r J ) is the time averaged intensity at r jt n and T are the running 
time and time delay, respectively, measured in units of the sampling interval r. 


In the experiments, the local intensity fluctuations were measured at two points 
only, distant by approximately one half the diameter of aperture H! . A spatial cor- 
relation index can then be defined as 


I<(\r x - r 2 1) = MaxC'(r 1 ,r 2 ,r), VT. (5.6) 

The correlation index runs from 0 to 1. It is high when the fluctuations at fq and 
r 2 are highly correlated. It is equal to one when the fluctuations at fq and r 2 differ 
only by a time delay, as one would expect if a transverse wave sweeps through the 
pattern. It drops to small values as the fluctuations at fq and r 2 become decorrelated. 
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Fresnel number F 


Figure 5.22: Spatial correlation index versus Fresnel number F. Comparison with 
Fig. 5.21 reveals the relationship between the spatial correlation and the number of 
defects. 
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The result is shown in figure 5.22. The interesting point is the correlation between 
the index I< and the number of defects (which measures the degree of spatial com- 
plexity), as the Fresnel number increases (i.e. , as the confinement of the system is 
relaxed). At first (Fig. 5.22), the correlation index drops, as one would expect. Near 
F ~ 4, where the first chaotic state appears, the correlation index drops drastically 
to a small value. In the range 4 < F <6, where the spatial complexity, as measured 
by the number of defects decreases (Fig. 5.21), the spatial correlation index is seen to 
increase (Fig. 5.22). The correlation index reaches a local maximum at exactly the 
same Fresnel number (~ 5.8) that leads to a local minimum in the number of defects. 
These observations seem to indicate that the loss of spatiotemporal coherence and the 
subsequent onset of turbulence in the phase-conjugate resonator is indeed mediated 
by the dynamics of the phase defects appearing in the wavefronts. 

5.6 Modes superposition 


Some of the transverse patterns observed are reminiscent of the common trans- 
verse electro-magnetic modes of a laser cavity. It is interesting to try to correlate our 
observations with the dynamical behavior of a few superposed cavity modes. The 
simpler motions observed at low Fresnel numbers could be interpreted as a sequential 
succession of a few modes of low order similar to the alternation described in refer- 
ence [99]. This interpretation seems to be confirmed by the high spatial correlation of 
these patterns. In addition, we have observed that the introduction of selective losses 
in the cavity can be used to suppress the oscillation of certain modes and lock the 
beam into an almost stationary pattern. At higher Fresnel numbers, when the motion 
becomes chaotic, the lack of spatial correlation seems to point to the coexistence of 
a number of uncorrelated modes as opposed to the alternation of modes of higher 
orders. The highly structured pattern observed at F ~ 5.8, with its higher spatial 
correlation and smaller number of defects, may however be due to the superposition 
or alternation of a small number of modes of higher orders. 

It is intriguing to ask whether a superposition of several modes with different 
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frequencies could produce patterns in which the defects would move around in a 
manner similar to the motions described in section 5.4. 

As an illustrative example we used a superposition of three Gauss-Hermite modes 
of the form 

E mn {x,y\t) = A mn U m ) exp }’ ^ 

where the Uj' s are Hermite polynomials, w is the beam “waist”, <f> mn = (m + n + l)fj 
and r) = tt/3 (arbitrary choice), 8 mn is the frequency detuning and A mn is the weight 

of mode mn [94]. 

We chose the three modes TEM 01, 10, and 20 with respective weights 1, 1, 1/3 
and respective detuning 1 / 2 , 1 , (\/5 + l)/2. A sequence showing the evolution of the 
pattern in time is shown in figure 5.23. The shaded regions are the bright spots and 
a and b indicate the position of the roots of Re(£)=lm(£)=0, where the defects 
are located. The arrows show the directions of their motions. Figure 5.24a shows the 
motion of the two conjugate defects projected on a plane. Starting at point A for 
defect a (solid dots) and B for defect b (empty squares), the following motions are 
observed. At first, defect a rotates slowly counterclockwise near the center the field, 
while defect b accelerates upward, disappears at the upper right, reappears at the 
upper left corner and swiftly moves down across the field. This motion is labeled 1 in 
Fig. 5.24a. Next (motion labeled 2), defect a , still moving slowly, changes direction 
and starts upward while defect b disappears at the bottom left and reappears at the 
bottom right. Defect a then accelerates out at the top left, reappears at the top right 
and rapidly crosses the field toward the bottom left, while defect b starts executing a 
slow clockwise rotation near the center of the field. In a certain sense, the two defects 
have exchanged their roles but with antisymmetrical motions. This type of motion 
continues, leading to the path diagram shown in figure 5.24b. 

Another interesting illustration given here is a scalar wave ip{r, t) combining the 
three modes TEM 00, 01 and 20 with fixed relative phases. For example, ip = 
cti{t)U 20 + MWoo + ia 3 {t)U 0 i, where ai(t)’s are real functions of time. This ex- 
ample shows two opposite charge defects moving apart along the x-axis as seen in 
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Figure 5.23. Sequence of patterns at equal intervals of time obtained by the superpo- 
sition of three Gauss-Hermite modes of order 01, 10 and 20 with respective weighting 
coefficients 1, 1, 1/3 and respective detunings 0.5, 1, (y/5 + l)/2. Letters a and b in- 
dicate the locations of the two conjugate defects and the arrows show their directions 
of motion. 
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Figure 5.24: Motions of the two conjugate defects seen in fig. 5.23 projected on a plane. 
The defects locations are solutions of Re£=lm£=0 where E is the field amplitude 
resulting from the superspoition of the three modes of Fig 5.23. Short time motion 
(see text for details) in top diagram, long-time defect trajectories in bottom diagram. 
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Figure 5.25: Computer-generated snapshots of an interferogram with a spherical ref- 
erence wave: before the nucleation of the defects (top diagram); while a defect crosses 
the vertex of the spherical wave (bottom diagram). 
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experiments with a Fresnel number F = 2.2 (figure 5.11). Figure 5.25 shows a com- 
puter generated interferogram of the wavefront with a spherical wave. Figure 5.25a 
shows the interferogram of the uniform wavefront with no defect. It displays typical 
Newton’s rings. A short while later, two defects appear at the center and move apart 
along the x-axis. Figure 5.25b shows the interferogram of the wavefront at the time 
when one of the defects crosses the reference spherical wave. The spiral fringes reveal 
the spiral nature of the wavefront. As the defect disappears at the boundary, the com- 
plete rings are recovered again. A similar behavior was observed in the experimental 
recordings for F = 2.2. Figure 5.26 shows snapshots of interferograms obtained by in- 
terfering the cavity beam with a spherical reference beam from the same laser. When 
one of the defects passes through the vertex of the spherical reference wave, it creates 

spiral fringes. 

These examples are presented for illustrative purposes only. They show that the 
superposition of a few low order modes, with appropriate weights and detunings, can 
produce a great variety of patterns and defect motions, some of them resembling 
qualitatively what is being seen in the experiments, but they do not prove that this 
is the correct representation. We will explore this modal superposition approach to 
represent the spatiotemporal dynamics in a PCR in the next chapter. 

5.7 Discussion and conclusion 

As mentioned in chapter 1, phase conjugation via four- wave mixing in BaTi0 3 is 
known to be itself unstable in some regions of parameter space, leading to a pulsating 
or chaotic phase-conjugate beam [22, 27]. These instabilities appear to be related 
to the existence and competition of several gratings in the crystal. In addition, self- 
phase-conjugation of pump two, which is readily observed in BaTi0 3 , produces a 
feedback to the laser which can also cause instabilities. Two different tests were 
made that lead us to believe that these two factors did not play a major role in the 
experiments. 

Four-wave mixing without the cavity (mirror M removed) and with an external 
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Figure 5.26: Experimental snapshots of the interferograms of wavefront containing a 
pair of defects with a spherical reference wave: before the nucleation of the defects (top 
diagram); while a defect crosses the vertex of the spherical reference wave (bottom 
diagram). 
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probe beam was tested for parameter values comparable to these used in the exper- 
iments. In all cases, the phase-conjugate beam was found to be quite stable. This 
appears to indicate that the observed dynamics is that of the phase-conjugate res- 
onator, which of course involves a possibly complicated nonlinear interaction with 
feedback between the resonator and the photorefractive crystal, rather than that of 
the four-wave mixing process. Furthermore, PC instabilities often require a rather 
high coupling constant (7/ ~ 6 in Ref [22]). In our crystal, the coupling constant was 
estimated to be between 2 and 3 only (see chapter 3). 

Tests were also made to detect the presence of self-pumping. With pump two alone 
incident on the crystal, strong self- phase-conjugation was observed. The presence of 
self- pumping is easily revealed by the usual beam fanning and internal reflections 
in the crystal. The feedback of the self- conjugated pump beam into the laser did 
cause some instability of its output. However, as soon as pump one was turned 
on, the fanning died out, the path of pump two straightened out inside the crystal 
and the self- conjugated beam intensity dropped to negligible values. Apparently, the 
presence of pump one did effectively erase all the extraneous gratings responsible for 
self-pumping. The same transverse patterns and dynamical behavior were observed, 
confirming that feedback due to self-pumping was negligible in our experiments and 
was not influential in our data. 

In summary we have studied some aspects of the spatiotemporal dynamics of a 
phase-conjugate resonator. The control parameter was the resonator Fresnel number. 
Time series of the local intensity fluctuations were recorded and analysed by means 
of spectral analysis and embedded phase space portraits. A number of occurrences 
characteristic of low dimensional dynamical systems were observed, i.e. in order of 
increasing control parameter: a bifurcation from a stationary to a periodic state, 
followed by more subharmonic bifurcations; the occurrence of a chaotic state char- 
acterized by the appearance of a continuous background in the power spectrum; the 
appearance of an additional frequency unrelated to the fundamental frequency of the 
periodic motion; the occurrence of intermittency; and the final development of what 
appears to be a fully turbulent state. 
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The occurrences of these events were found to be quite reproducible, providing 
some evidence for the presence of a deterministic chaotic state of the cavity field. 
Further evidence is given by calculating the correlation exponent and the Kolmogorov 
entropy at a Fresnel number leading to a chaotic time series. 

The spatial complexity of the pattern was studied by video recording and by 
interferometry, revealing the leading role played by the wavefront dislocations or vor- 
tices. The number of phase defects and the spatial correlation index of the wavefront 
were observed to follow similar trends as the transverse confinement of the system 
was relaxed. This suggests that the spatiotemporal chaos of the cavity field may be 
mediated by the phase defects. 
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Chapter 6 


A model of spatiotemporal 
dynamics in PCRs 


6.1 Introduction 

In this chapter we present a truncated modal expansion approach to model the 
spatiotemporal dynamics of a PCR. Nonlinear dynamical problems that include the 
transverse dimensions and are described by systems of partial differential equations 
(PDEs) remain a theoretical challenge and do not benefit from well established 
experimental methodologies [100] as do systems of ordinary differential equations 
(ODE) [60, 101]. Furthermore, their numerical analyses require large computational 
budgets [56, 102]. We will test a relatively simple mathematical scheme to represent 
the transverse dynamics of the PCR without requiring prohibitive computational 
budgets. 

This chapter has two aims. The first is to describe and establish the soundness 
of a somewhat arbitrary yet sensible modal decomposition approach to study the dy- 
namics of PCRs with modest Fresnel numbers. This is done by comparing numerical 
results with actual experimental data. The hope is that this comparison, although 
complicated by the difficulty in measuring experimental control parameters exactly 
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and by the problematic effect of stochastic noise, will prove sufficiently convincing to 
confirm the usefulness of the approach. The second is to illustrate the relationship 
existing between the dynamics of the optical vortices observed in PCRs [87, 88, 86] 
and the spatial correlation of the oscillating beam. 

The chapter is organized as follows: The PCR’s equations and boundary con- 
ditions within approximations compatible with experimental conditions are stated 
in section 6.2.1. Arguments justifying the modal expansion approach are given in 
section 6.2.2 and the method is then used to arrive at a set of modal amplitude equa- 
tions. In section 6.3, a particular set of parameter values is chosen and the system 
is integrated, using the off-Bragg detuning as a control parameter. Different dynam- 
ical behaviors are identified from local time series, power spectra and phase space 
portraits. Examples of vortices’ motions and spatial correlation maps are given in 
section 6.3. Experimental results are presented and discussed in section 6.4 and a 
summary is given in section 6.5. 


6.2 The modal expansion model 


6.2.1 Phase-conjugate resonator 


The PCR is sketched in fig. 6.1. Its active element is an externally pumped 
photorefractive medium acting as a PCM. Maxwell’s equations for the optical field in 
the medium lead to the following coupled equations (see Appendix A): 
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Figure 6.1: Four-wave mixing in a photorefractive phase-conjugate resonator (PCR): 
A u A 2 are two gaussian pump beams; A 3 , A 4 are the left-going and the right-going 
beams inside the PCR; mirror M allows feedback of field back into the phase-conjugate 
mirror (PCM); apertures 1R ,H 2 confine the transverse dimensions of the propagating 

beams. 
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where, in spirit of the slowly varying envelope approximation, the field in the medium 
is written as 

4 

E(x,y,z,t) = ^Ai(x,y,z,t)e i(ui - k '^ + ex., (6.5) 

1=1 

v — cjn is the speed of light in the medium, n is the average refractive index, kj are 
unit vectors pointing in the direction of the four wavevectors in the medium. The 
wavenumber is kj — u>/v, where u> is the angular frequency of the pump and a is the 
linear absorption coefficient. Vx is the transverse gradient operator. The parameter 
b is the wavevector mismatch (momentum mismatch along the z-axis) measuring 
the departure from the Bragg condition. This off-Bragg parameter is introduced as 
a means of varying the amount of phase transfer between the waves mixed in the 
nonlinear medium and is later taken as the control parameter. 


The equation for the phase grating Q is derived from the standard band trans- 
port model of Kukhtarev [103] and reads, in the single grating approximation (see 
Appendix B), 

9Q(x,y, z,t) 


dt 


+ CQ(x,y,z,t) = 


( 6 . 6 ) 


(j ( Ai(x,y, z, 0 ^ 4 1 3 ') y_ ) z t d~ ■A 2 ( z, y, tyAs(x^ y, £, t ) 

V l Q (x,y,z,t ) 


where C = [E 0 + i(E q + E D )\/[E a + i(E M + E D )\. 


The boundary conditions for the wave equation in the medium are determined by 
the cavity geometry. The transverse confinement of the cavity is controlled by two 
square apertures Hi and H 2 of sides 2 a pcm and 2a m respectively. The first aperture is 
located near the phase- conjugate mirror in the front focal plane of the lens (fig. 6.1). 
The second one is in the back focal plane of the lens next to the planar mirror. The 
transverse confinement is measured by the cavity Fresnel number defined as 


^ 4a m a pcm / A/, (6.7) 

where A is the wavelength in vacuum and / the focal length of the lens. 

One cavity round trip, starting at the PCM, consists of a truncation by H t and 
imaging in a 4-f afocal system with a pupil H 2 . Thus, the amplitude A 4 is proportional 
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to the truncated convolution of A 3 with a point-spread-function which is the fourier 
transform of the aperture H 2 . According to Fourier optics [104], the field at the plane 
of the mirror M, £(x, y, -I + , t), is proportional to the convolution of the Fourier 
transform of the field A 3 (x, y,0 + , f) at the exit plane of the crystal and the Fourier 
transform of the square aperture Hx, i.e. 


L e'koL 
£(x,y,-L + ,t+-) - -jj 


F{A 3 {x , y, 0 + , t)} * F{ rect(^— )rect( ; V 


2 a 


pcm 


2a 


■» 


, ( 6 . 8 ) 


pcm 


where rect(x/2) =1 for — 1 < x < 1 and = 0 otherwise. The field is then mirror 
reflected and truncated by the aperture H 2 such that 


£(x,y,-L + ,t + -) = r£(-x,-y,-L + ,t + -) • rect( 
C C 


X )rect(-f ), 


2a 


< 2 a r 


(6.9) 


where r is mirror M’s reflectance coefficient. 

The field £{x,y, is a multiplicative of the inverse Fourier transform of 

£{x,y , -L + ,t + \ ), i.e. 


, _ 2 £, 

£{x,y, 0 ,< + — ) 


e TlF-'{£(x,y,-L+ 1 t + -)) ( 6 - 10 ) 

JAJ C 

e ' 2k ° L -A 3 (-x,-y,0\t) ■ r ect(^)rect(^) 


A 2 / 2 


2a 


pcm 


2 a 


pcm J 


* 4a^sinc(2a m x)sinc(2a m i/)] , 


where sinc(x) = sin(7rx)/(7rx). 

Thus, we can write the input field A 4 (x, y, 0 + , t) in terms of the outgoing phase 
conjugate beam A 3 (x, y, 0 + , t) in the following way: 


9 J , p i2 k 0 L - f-\-a pc m 
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Using the resonator Fresnel number F = 4a m a pcm / Xf and normalizing the trans- 
verse coordinates to the aperture size a pcm , we find 

*2d L j 

A 4 (x,y,0,t + —) = ~r(F/2j 2 e ,2k ° L JJ ^ A>(x', y\ 0, <)sinc [F{x' + x)/2] (6.12) 

sine [F(y' + y)/ 2] dx'dy', 

where x — x/a pcm ,y = y/u pC m ( x and y are the real physical dimensions). 

6.2.2 Modal decomposition 

The direct integration of equation 6.1- 6.6 requires a large computational bud- 
get. Furthermore, calculations taking only one transverse dimension into account, 
although showing interesting dynamical behaviors (e.g. ref. [105]), fail to reveal the 
crucial role played by the phase defects in mediating the spatiotemporal dynamics of 
the PCR [87, 88, 86], 

Modal decomposition is a powerful tool for studying systems of PDEs [106]. This 
approach allows one to reduce an infinite dimensional problem to a finite (hopefully 
low) dimensional problem described by a system of ODEs for the modal coefficients. 
The success of the method relies on the assumption that the dynamics is fairly well 
described by the interaction of a small number of active spatial modes. The main 
difficulty resides in the choice of appropriate basis functions for the expansion. They 
must accurately represent the spatial structure of the active modes. However, this 
structure is generally difficult to identify a priori and one often relies on the choice 
of a “sensible” basis rather than trying to identify the exact optimal basis. 

In optical resonators with relatively small Fresnel numbers, this assumption ap- 
pears to be justified. The dynamics of these systems is indeed dominated by a few 
active modes and the Gauss- Laguerre or Gauss-Hermite modes of the empty res- 
onator seem to be reasonable basis functions to replace the optimal ones [75], Of 
course such a practice introduces a degree of arbitrariness since a small number of 
modes of the sensible basis may not be sufficient to adequately represent the spatial 
structure. Nevertheless, the resulting enormous reduction of computational budget 
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and the physical insight that can be gained justify further tests of the validity of this 

approach. 

The cavity fields A 3 and A 4 in the nonlinear medium (0 < z < d) are thus 
expanded in series as 


A 3 (x,y,z,t ) = *;*<>) ( 6 - 13 ) 

771,71 

A 4 (x, y, z,t) = '%2 hmn{z, t)u m {x , z; z 0 )u n {y , z; z Q ) (6.14) 

771,71 

and the Gauss-Hermite modes are chosen as a sensible set for the decomposition: 


(•£ ? Z , Zq ) 



exp 
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u; 2 (z) 
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■f i(m + 1/2) tan l (zlz 0 ) , 


( 6 . 15 ) 


where H m are Hermite polynomials. The beam radius is w{z) = uj 0 (1 + z 2 /z 2 ) 1/2 , its 
radius of curvature is R{z) = z(l + z 2 /z 2 ). z„ = Tnv 2 0 n/X is the Rayleigh distance 
an( j _ 27 rh/A is the wavenumber in the medium. The waist w 0 is located at the 

entrance face of the PCM (z = 0). 


The pump beams are Gaussian, with a waist w p chosen appreciably larger than 
w 0 , and for simplicity are assumed to make a very small angle with the z-axis. 1 hus, 
the pump beams are 


A 1 {x,y,z,t) = /i(-M)t/(x,z)£%,z) ( 6 - 16 ) 

A 2 (x, y, z, t) = / 2 (z, t)U*{x, z)U*(y , z) (6.17) 


where U (x, z) = u 0 (x, z; z p ) and z p — nw p n/ A. 

The series in equations 6.13, 6.14 are infinite and the problem remains infinite 
dimensional. Because of the cavity transverse confinement however, we expect that 
only the modes with indices smaller than some upper bound will take a significant 
part in the dynamics. The cavity Fresnel number defined by equation 6.7 is a measure 
of the amount of transverse spatial information that the cavity can accommodate ( F 
is the space-bandwidth product of the cavity). It is thus reasonable to assume that 
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only the modes with indices smaller than F will play a dominant role in the wave 
interaction in the PCM since they are the modes which are most likely to survive in 
the cavity. It must be stressed however that choosing which modes to include in the 
dynamics is somewhat arbitrary. Only a comparison with experiment can provide 
an a posteriori justification of this choice. The important point is that the series of 
equations 6.13, 6.14 can often be limited to a small number of terms while still giving 
results which agree at least qualitatively with observations, as will be shown in a later 
section. 


Equations 6.13- 6.17 are substituted in 6.1- 6.6 and use is made of the fact that 
the Gauss-Hermite modes are solutions of the Helmholtz equation: 

(iv 2 T d \ 

1 - 2 *- + Yz ) Um ^ x ' z ' Un ( y ' z ) = °> ( 6 . 18 ) 

and satisfy the biorthogonality relation [94] 
r y+ oo 

J dxdyu m (x, z)u n {y , z)u p (x , z)u*[y, z) — S mp S nq . (6.19) 

A final set of differential equations is then obtained for the modal amplitudes: 
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( 6 . 21 ) 

( 6 . 22 ) 

(6.23) 


where 


hmn{z,t)=ff dxdyQ{x,y,z,t)u m {x,z)u n (y,z)U’(x,z)U*(y,z) (6.24) 

J J — OO 


is the overlap integral of the pump and the projection of the grating Q onto the mode 
(m, n). 
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The grating equation 6.6 is likewise projected onto the cavity modes to obtain 

(6.25) 


r a / w ( M ) + chmn ( z , t) = 


dt 


^ 0 C{e xb 2 f 1 {z,t)Y,fw(^i)Gklrnn{z,t) + f 2 (^i)'Lf 3 klGklrnn{z,t)} : 
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, x _ f f + °° u k{ x > ^) u r(y^) u "i( x > z )u n (y, z)\U(x, z)| \U(y,z)\ , , 26) 

Gklmn{z, t) — J J ^ I 0 (x,y,Z,t) 

and the total intensity I 0 (x,y,z,t) is calculated using equations 6.13, 6.14, 6.16, 6.17. 

For the case of strong pumps and large pump width (we can call this the zeroth 

order approximation), G klm n in equation 6.26 becomes approximately a delta function 

8 km Si n such that the grating equation becomes 

r d/w(z,j ) +C h mn (z,t) = lo C (6-27) 

at 

{e lilZ /l(z, t)/4 mn (z, ^) + f2{ z ,i)f3mn }i 


It seems that as the coupling terms in the sums of the right hand side of equation 6.25 
have been taken out in this approximation, the projection of the grating Q onto each 
mode (m,n) contributes mainly to the corresponding mode (m,n). However, cou- 
pling still occurs through the factors f\ and /£ which are linked to all participating 
modes as shown in equations 6.20 and 6.21. Thus, in spite of the fact that coupling 
through G klmn is eliminated in the strong pumps and large pump width approxima- 
tion, coupling between modes remain in equation 6.27 due to factors f\ and / 2 , unless 
if the pumps are so strong that f\ and 3-re about constant. We may regard this 
so-called zeroth order approximation as appropriate to the spatiotemporal dynamics 
with weak coupling in PCR. 

The integration of equations 6.20- 6.26 is greatly simplified if the interaction length 
in the PCM is shorter than all the beams’ Rayleigh distances: 

d<z 0 ,2p- (6.28) 


In this case, diffraction can be neglected in the medium, and the interacting beams 
modes, normalized to the aperture size a pcm , take the simple form: 
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with /? 0 = y/2a vcm /w 0 . The pumps are gaussian with a parameter /3 p = \/2a pcm /w p . 

The boundary conditions for the modal amplitudes at 2 = 0 are found from 
equation 6.12 using the modal decompositions 6.13, 6.14. They read, without the 
unimportant phase factors in eq. 6.12, 


/4mn(0,0 = rF 2 '£f 3kl (0,t)I km I ln , (6.30) 
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TC k (x) = j dx'u k (x')smc[F(x' + x)/2]. 


(6.31) 

(6.32) 


The algorithm of the adiabatic elimination method described in section 3.1.2 was used 
to solve the system of equations 6.20- 6.26 with the boundary conditions 6.30- 6.32. 
Some results are illustrated in the next section. 


6.3 Numerical results 


Phase transfer between interacting beams in the four- wave mixing process appears 
to be a desirable condition for the observation of non-trivial dynamics with reason- 
able values of the coupling parameter. For example, the onset of instabilities in ring 
resonators with injected signals has been shown to depend sensitively on this param- 
eter [24]. In a material having a real coupling parameter and no external field, this 
phase transfer can be the result of an angular mismatch of the pump beams. For 
the simulation presented in this section, the control parameter was chosen to be the 
phase mismatch bd (off- Bragg parameter), where b is the momentum mismatch along 
z and d is the length of the nonlinear medium. 

We have shown in the last chapter that the dynamics of the PCR strongly depends 
on the cavity Fresnel number. Increasing the Fresnel number relaxes the transverse 
confinement of the system and allows higher order modes to take part in the dynamics. 
As a result, the spatial complexity of the beam increases. In particular, the number 
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of vortices increases and their motion becomes more complex, eventually leading to 
chaotic states and, in the limit of large Fresnel numbers, to spatiotemporal turbulence. 
The example shown in the next section is for a Fresnel number F = 3.5. This is small 
enough so that only a few modes are expected to take part in the dynamics, yet it is 
large enough to exhibit a wide range of behaviors. 

Experiments carried out with Fresnel numbers between 3 and 4 reveal a transverse 
pattern in which two pairs of defects of opposite charges repetitively nucleate, move 
along symmetrical trajectories, spending much time in some circular zone centered on 
the beam and disappear at boundaries or collide and annihilate each other [88]. Based 
on these observations and in order to keep the computational budget to a minimum, 
it was decided to limit the expansions in eqs. 6.13, 6.14 to four modes only, namely, 
(m,n) = (0,0), (1,1), (2,0), and (0,2) with a beam waist w Q = a pcm = 1. 

As already mentioned, there is a certain degree of arbitrariness in this choice. How- 
ever, with a Fresnel number of 3.5, it is unlikely that modes with indices higher than 
3 can survive in the cavity and play a significant role in the dynamics. Indeed, when a 
mode of higher index (8, 0) was initially added to the series of (0, 0), (1, 1), (2, 0), (0, 2), 
the simulation showed that the amplitude of mode (8, 0) decays and vanishes with 
the other transients. Another reason for limiting the expansion to these particu- 
lar four modes is that it allows the calculation of the vortices’ positions to be done 
algebraically rather than numerically, thus saving CPU time. Finally, simulations 
run with different choices of modes, although giving results which differ in details, 
revealed the same generic range of behaviors. 

The other parameters of the model are chosen as follows: 7 d = —10, ad = 

0.15, & = y/2,p p = V2/10 ,E d = 1 kV/cm,E q = 5kV/cm,E M = mkV/cm,E 0 = 0, 
and the pump beam intensity coefficients are f?(z = 0) = j\{z = 0) = 0.5. As an 
initial condition, it is assumed that the four modes are excited from white noise and 
have initial modal amplitudes hmn{t = 0) = 10- 4 . 
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6.3.1 Local intensity fluctuations 

The local intensity of the cavity field A 4 measured at the exit face of the PCM was 
chosen as a variable because it is easily accessible experimentally. The local intensity 
h(x 0 ,y 0 ,d) was calculated at (x 0 ,y 0 ) = (0.75,0.15) in a region where experiments 
indicate that the intensity fluctuations have a deep modulation. Figure 6.2 shows 
short sections of time series of the local intensity fluctuations for different values of 
the off-Bragg parameter. Time in fig. 6.2 is normalized to the relaxation time r' of 
the photorefractive grating (r' = \[E 0 + i(E M + E D )]/[E 0 + i(E q + E d )}\t„ see eq. 6.6). 

After the transients have vanished and a stationary state has been reached, the 
local extrema of these time series were recorded and used to construct the bifurcation 
diagram of fig. 6.3. For small Bragg mismatch (\bd\ < 2.4) the motion is periodic. 
Unstable oscillations occur in a wide range of parameter values (2.4 < | bd | < 3.6). 
In most of this range the motion appears to be quasiperiodic. Some regions of the 
bifurcation diagram seem to be more chaotic than others and there exist narrow 
windows in which the motion appears to be periodic, usually with a large number 
of subharmonics, for \bd\ > 3.6, the motion is periodic again and eventually, with a 
sufficiently large mismatch, the only stable state is I 4 = 0. 

The time series of fig. 6.2 were chosen to illustrate these various types of motion. 
The corresponding power spectra shown in fig. 6.4 give a fair indication of the dy- 
namics of the system. Additional information can be gained by reconstructing pseudo 
phase space portraits such as these shown in fig. 6.5, which are plots of the imaginary 
part of the field A 4 (x 0 ,y 0 ,d) versus its real part. 

At bd = —2.3, the motion is clearly periodic with a frequency fj slightly larger 
than 1/2 of the relaxation rate of the grating (fig. 6.4a), and the phase space portrait 
(fig. 6.5a) is a limit cycle. At bd = —2.45, a second characteristic frequency f 2 , about 
1/10 of the grating relaxation rate, appears and modulates the time series (figs 6.2b 
and 6.4b). For this value of the off-Bragg parameter, the two frequencies are (rela- 
tively ) incommensurate and the motion is quasiperiodic with the phase space portrait 
shown in fig. 6.5b. At bd = —2.55, the low frequency modulation deepens (fig. 6.4c) 
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Fieure 6.2: Time series of local intensity J 4 for different values of the off-Bragg pa- 
rameter bd showing a sequence of dynamical states including periodic, quasipeno ic 

and chaotic states. 
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Figure 6.3: Bifurcation diagram of the local intensity I 4 at a point of the beam cross- 
section in the exit plane (z = d) as a function of the off-Bragg parameter bd. Other 
parameter values used are: ~fd = -10, ad = 0.15, F = 3.5 ,/3 0 = 1 ,/3 p = 0.1, /^(z = 
0) = fi(z = d) = 0.5, = 0) = IQ" 4 , (m, n) = (0, 0), (1,1), (2, 0), (0, 2). 
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Figure 6.4: Power spectra of the time series of fig. 6.2 illustrating the dynamics as 
the off-Bragg parameter bd is varied, (a) a periodic motion with a single frequency, 
(b) quasiperiodic motion with two (relatively) incommensurate frequencies, (c) same 
with more energy in the side bands, reflecting the deeper modulation of fig. 6.2c, (d) 
broadening of the two frequency peaks and their harmonics with the appearance of 
an additional low frequency leading to chaotic motion, (e) two frequency peaks of 
commensurate ratio corresponding to frequency locking at period 13. 
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Figure 6.5. Phase space portraits corresponding to the local time series of fig. 6.2 with 
the real part of the field amplitude A 4 plotted versus its imaginary part. It shows 
(a) a limit cycle, (b) a two frequency torus, (c) a thicker two frequency torus, (d) a 
diffuse broken torus, (e)a 13 loops limit cycle indicating frequency locking. 
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and more energy shifts toward higher harmonics of U and linear combinations of the 
two main frequencies (fig. 6.4c). The phase space portrait still shows some structures 
although it will eventually densely fill a region of phase space since the two frequencies 
are still incommensurate. At bd = -2.75, the time series is more irregular, possibly 
chaotic (fig. 6. 2d). There are still two main frequencies but the spectrum is broadened 
and shows a number of additional peaks indicating the presence of a third frequency 
about 1/100 grating relaxation rate (fig. 6.4d). The phase space portrait of fig. 6.5d 
is diffuse and shows no apparent structure. At bd = -3.18 a phenomenon akin to 
frequency locking occurs. The two larger frequencies pull each other so as to become 
rational. The very low frequency, which seems to have appeared only to allow this 
locking to occur, disappears. The spectrum of fig. 6.4e shows a ratio fi/f 2 =13. The 
periodic motion with 13 subharmonics is confirmed by the closed loop phase space 

trajectory of fig. 6.5e. 

6.3.2 Vortices and spatial correlation 

The local intensity fluctuations offer a convenient means of comparing model pre- 
dictions with experimental data but this information is not sufficient to fully charac- 
terize the spatiotemporal dynamics. The aim of this section is to illustrate the role 
played by the phase defects (vortices) in mediating the dynamics. 

The vortices appear where the real and imaginary parts of the field amplitude 
vanish simultaneously. For example, fig. 6.6 shows a snapshot of the beam ampli- 
tude cross-section |A 4 (x,y,d)| for bd = -2.45, which exhibits four dark spots. The 
corresponding phase contour diagram of fig. 6.7 clearly identifies these spots as two 
pairs of vortices of opposite charges (±1). In time, these vortices nucleate, move 
around, annihilate each other or disappear at boundaries. This is the motion of the 
vortices which gives rise to the local intensity fluctuations described in the previous 
section. This section gives a more detailed illustration of the spatial aspect of the dy- 
namics and attempts to identify the spatial features responsible for the characteristic 
frequencies observed in the local times series. 
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Figure 6.6: The profile of the absolute value of the field amplitude 
noles in the beam cross-section. 
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Figure 6.8: Time evolution of the modulus of the field amplitude A 4 along a line x = 1 
for different values of the off-Bragg parameter bd , showing (a) a synchronized periodic 
motion, (b) a periodic motion with transverse modulation, (c) a quasiperiodic motion, 
(d) a turbulent flow, (e) another periodic flow. 
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Figure 6.8a to e shows the temporal evolution of the modulus of the field amplitude 
A a { X = 0.75, y, d ) along a line across the beam for the same off-Bragg parameter values 
used in figs. 6.2, 6.4, 6.5. For a small offset ( bd = —2.3), all points across the line 
execute synchronized periodic motions with a period a bit smaller than two grating 
time constants (fig 6.8a). In 2-D, the phase of the wavefront is observed to breathe 
periodically together with the amplitude but the phase gradients produced are not 
steep enough to trigger the nucleation of defects. 

At bd = -2.45 (fig 6.8b), the periodic motion is transversely modulated with a 
period of the order of 10 grating time constants. This spatial motion gives rise to 
the quasiperiodic time series of fig. 6.2b. In 2-D, the phase gradients are locally 
steep enough to tear the wavefront at locations where two pairs of defects of opposite 
charges nucleate. The four vortices then travel across the wavefront, and pairs of 
opposite charge collide and annihilate. The higher frequency of the time series seems 
to correspond to the recovery rate of the local wavefront after a vortex has moved 
through it, while the lower frequency corresponds to the full cycle of vortex nucleation, 
motion and annihilation. At bd = -2.55, the motion is similar to that just described 
but with deeper modulation (fig. 6.8c, which is plotted with a different time scale). 

bd = —2.75 (fig. 6.8d), the transverse modulation is irregular, possibly due to 
the occurrence of a third incommensurate frequency. In 2-D, in contrast with the 
previous case where the vortices’ trajectories were confined to a fairly well defined 
area, they now seem to visit the whole beam cross-section irregularly. At bd = -3.18 
(fig. 6.8e) the motion is periodic again. The transverse modulation has a period of 
~ 23 grating time constants, which is 13 times the period of ~ 9/5 time constant of 
the fast oscillations measured in fig. 6.8a. 

The remainder of this section gives a more detailed illustration of the relationship 
between the vortices’ trajectories and the spatiotemporal coherence of the beam. 
The dynamics of the fluctuating beam can be characterized by a correlation index 
distribution calculated from the spatial correlation function defined in eq. 5.4. 

To illustrate the usefulness of this concept, the correlation index distribution with 
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respect to the origin, <7(0,0; x,y), was calculated for the five off-Bragg parameter 
values used in the previous examples. The spatial correlation diagram thus obtained 
can then be compared with the vortices’ trajectories. For this, the positions of the 
vortices were calculated at regular intervals of time by solving for R eA 4 (x,y,t) = 
lmA 4 (x,y,t) =0, and they are represented by dots on a 2-D plot. Note in passing 
that, with the particular choice of modes adopted for the expansion of A 4 , these 
positions could be solved for algebraically. 

Figure 6.9 shows the correlation index distribution for a small value of the off- 
Bragg parameter (bd = -2.3). No vortices appear in the field but the amplitude 
and phase of the wavefront oscillate periodically as if waves were travelling along an 
annular area centered on the beam. 

At bd — —2.45, the motion is characterized by the nucleation of two pairs of 
defects at two diagonal ends. The vortices move along circular trajectories and the 
members with opposite charges from each different pair collide and annihilate near 
the two diagonal ends on the other side. This motion is then repeated periodically 
with alternating directions. A sequence of interferogram snapshots illustrating this 
motion is shown in fig. 6.10. A plot showing the distribution of vortice’s positions in 
time is shown in fig. 6.11a. The vortices’ trajectories are clearly confined to a narrow 
annular area in the beam. Figure 6.11b shows the corresponding correlation index 
distribution. Its particular shape can be explained by the fact that the intensity fluc- 
tuations at some location in the beam are due to the passage of a vortex nearby, that 
no vortex ever appears near the origin (the fixed point for calculating C) and that the 
vortices’ trajectories are confined to a narrow annular region of space. It is only near 
this region that the correlation index is expected to drop sharply. The correlation in- 
dex distribution was also calculated with a fixed point chosen at a location frequently 
visited by a vortex. The relationship between the distributions of fig. 6.11a and b 
is analogous to that between a photographic print and its negative. The correlation 
index is high in the annular area visited by the vortices and low everywhere else. This 
confirms that the four vortices appearing in the pattern are highly correlated, as is 
expected since they nucleate simultaneously and produce intensity fluctuations near 
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bd= -2.3 
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Figure 6.9: Correlation index distribution calculated with respect to the center point 
shows that the spatial coherence remains high in most parts of the wavefront except 
near an annular area where phase instability is found to occur. 
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Fieure 611: (a) Trajectories of defects in the wavefront for the off-Bragg parameter 
U = -2.45, localized on an annular ring, (b) The correlate mdex with respect to 
a point at beam center drops sharply near the defects’ trajectories. 
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their trajectories which are only delayed in time. A qualitatively similar behavior of 
the vortices was observed for bd = -2.55, as illustrated in fig. 6.12a and b. 

An apparently chaotic time series was observed for bd = -2.75. Figure 6.13a shows 
that the vortices’ position distribution is diffuse, indicating that their trajectories do 
not remain confined to narrow regions of space. In fact, the defects can visit about any 
place in the beam except near the origin and this is only because of the peculiar set of 
modes chosen for the decomposition. The corresponding correlation index distribution 
shown in fig. 6.13b is low in all regions visited by the vortices. Figure 6.13c shows the 
correlation index distribution with a fixed point = (1,0.3) frequently visited 

by a vortex. It shows that even though the vortices trajectories have become more 
random, the intensity fluctuations in the regions they visit remain highly correlated. 

At bd = -3.18, an apparent frequency locking was observed in the local intensity 
time series, leading to a quieter dynamics. For this value of the off-Bragg parameter, 
the wavefront phase and amplitude oscillate periodically but no defect nucleates. The 
regions in space where these fluctuations occur are identified as two concentric rings 

in the interferogram snapshots of fig. 6.14 and in the corresponding correlation index 
distribution of fig. 6.15. 


6.4 Experimental results 

The experimental apparatus has been described in detail in section 5.1, The 
PCM is an externally pumped single crystal of BaTi0 3 . The source for the pumps 
is a single mode Ar ion laser (A = 514nm) optically isolated from the setup by a 
Faraday rotator. The cavity ends with a planar dielectric mirror (R=0.95%) and 
contains a 16cm focal length lens. Two intra-cavity circular apertures of diameter 
d P cm = 0.79mm and d m =0.38mm give a Fresnel number F = 3.7. The off-Bragg 
parameter was varied by tilting the mirror directing the pump A l toward the PCM 
and its values were calculated as the product of the momentum mismatch along the 
cavity axis with the interaction length. The change in pump overlap due to the tilting 
was measured and found negligible. 
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Figure 6.13: (a) At bd = —2.75, the defect trajectories become more diffuse and 
occupy larger regions of space in the wavefront, (b) The correlation index distribution 
with respect to the center point drops to small values at places where the defects 
concentrate, (c) The correlation index distribution with respect to a point (1,0.5) 
shows a high spatial coherence in all areas visited by the defects. 
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Figure 6.15: The correlation index distribution with respect to the center point show 
a structure with two concentric rings. The regions of low values of the correlation 
index are the places where fringe bendings are observed. 
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The aim of the experiments is to show, by way of illustration, that the range of 
behaviors predicted by the model and described in section 6.3 is indeed observed in 
the physical system. It would be unrealistic to expect a quantitative match between 
the experimental data and the results of numerical simulations if only because it is 
not possible to control experimentally which modes are actually taking part in the 
dynamics. A qualitative match, however, would at least indicate that the truncated 
modal expansion approach gives a fair representation of the actual dynamics of the 

PCR. 

Figure 6.16 shows an interferogram snapshot of the beam exiting the PCM. It 
exhibits four vortices located at the corners of a quadrangle in an arrangement similar 
to that predicted by the model in figs. 6.6 and 6.10. 

The off-Bragg parameter was scanned between two extreme values at which the 
momentum mismatch is large enough to reduce the gain below cavity losses and pre- 
vent oscillation. Local time series were recorded within that range and analyzed using 
power spectra and delayed-time phase space portraits. Near the center of the cov- 
ered range of parameter values, a region where the motion was clearly periodic could 
be identified. On both sides of this region, similar sequences of various dynamical 
behaviors were observed, up to large mismatches at which the cavity stops oscillat- 
ing. Before reaching these limits, simple periodic oscillations at reduced amplitudes 
were again observed. The sequence of behaviors just described is exactly what the 
bifurcation diagram of fig. 6.3 would predict (note that within the approximations of 
the model, this bifurcation diagram is symmetrical around bd = 0, but in the actual 
setup this symmetry is broken). By analogy with the model, the origin of parameter 
space (bd r = 0) was chosen as the midpoint of the central region showing periodic 
oscillations. Other values of the off-Bragg parameter mentioned below are relative to 
this point. Some examples of time series are shown in fig. 6.17 and the corresponding 
power spectra are displayed in fig. 6.18. These examples were chosen to illustrate the 
variety of dynamical behaviors exhibited by the PCR and to compare them with the 
characteristic behaviors predicted by the model. 

Figures 6.17a, 6.18a show the periodic motion expected for the off-Bragg parameter 
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Figure 6 . 16 . Interferogram snapshot at the output of the PCM revealing the presence 
of two pairs of vortices of opposite charges (experimental result). 
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Figure 6.17: Local time series of I 4 measured at a point in the beam cross-section 
showing the variation of the dynamics for different values of the off-Bragg parameter: 
(a) simple periodic motion, (b) irregular motion, (c) frequency locking, (d) quasiperi- 
odic motion, (e) another quasiperiodic motion. 
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Figure 6.18: Power spectra of the local time series of fig. 6.17 showing (a) a funda- 
mental frequency for a simple periodic motion, (b) three frequencies and a continuum 
for an irregular motion, (c) a frequency and it subharmonics for the frequency locked 
motion, (d) two incommensurate frequencies for a quasiperiodic motion, (e) another 
quasiperiodic motion with at least two incommensurate frequencies in the spectrum. 
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bd = 0. The corresponding delayed-time phase space portrait (fig. 6.19a) shows 

broadening which may give an indication of the level of stochastic noise. 

Figures 6.17b, 6.18b obtained with bd T = 0.65 show a chaotic state with an ir- 
regular time series and a spectrum containing at least two or three incommensurate 
frequencies, including the main frequency of about 0.14Hz, and a broad continuum, 
a signature of chaos. The exponential decay of the power spectrum (linear on the 
semilog scale of fig. 6.18) may be indicative of deterministic chaos. The corresponding 
phase space portrait of fig. 6.19b is diffuse and essentially featureless. 

Figures 6.17c, 6.18c give a typical example of frequency locking obtained at bd T = 
1.29. The locking at period three is confirmed by the closed-loop phase space trajec- 
tory shown in fig. 6.19c. Note that the width of this limit cycle is about the same as 
that of fig. 6.19a and is probably due to stochastic noise. 

At bd r = 1.94 (figs. 6. 1 7d, 6.18d) the motion is quasiperiodic with a spectrum 
containing only two incommensurate frequencies and their linear combinations. An- 
other quasiperiodic motion obtained at bd r = 2.16 is shown in fig. 6.17e, 6.18e. Here 
the lower frequency appears as sidebands of the higher one. For larger values of the 
off-Bragg parameter, e.g. bd T > 2.37, the motion is periodic again with the same 
dominant frequency of ~ 0.14Hz but with decreasing amplitude. When bd r > 3, the 
gain is lower than the cavity losses and no oscillation occurs. 

6.5 Summary and conclusions 

A model based on a truncated modal expansion of the cavity modes has been 
tested to describe the spatiotemporal dynamics of a phase- conjugate resonator. The 
Gauss-Hermite modes of the empty cavity are chosen as sensible basis functions to 
represent the fields in the nonlinear medium. Although the choice of the modes 
which take part in the dynamics is arbitrary in this method, justifications are given 
for truncating the series at a mode index of the order of the cavity Fresnel num- 
ber. Verification of the validity of this choice are later given by comparison with 
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Figure 6.19: Time delayed phase-space portraits obtained from the data of the time 
senes in fig. 6.17 showing (a) a limit cycle, (b) an irregular distribution, (c) a cycle- 
three motion, (d) smeared trajectories for the quasiperiodic motion 
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experimental data. Numerical examples, using the olt-Bragg parameter (momentum 
mismatch) as a control parameter, reveal a rich variety of dynamical behaviors illus- 
trated by the local intensity fluctuations and ranging from simple periodic oscillations 
to quasiperiodic motions and chaotic states. 

The modal decomposition method is particularly well suited to studying the spa- 
tial aspect of the dynamics without requiring prohibitive computational budgets. 
Numerical results indicate that the spatiotemporal dynamics of the PCR is mediated 
by the nucleation of pairs of defects of opposite charges in the beam and by their sub- 
sequent motion and annihilation. Maps of the spatial coherence function of the beam 
are found to be templates of the corresponding maps of the vortices’ trajectories, 
establishing a strong correlation between these two quantities. The loss of temporal 
coherence and the onset of temporal chaos in the local intensity fluctuations is like- 
wise correlated to the loss of spatial confinement of the vortices’ trajectories and to 
the loss of spatial coherence. 

Experimental data obtained with a PCR, using an externally pumped BaTiOs 
photorefractive PCR, reveal a range of dynamical behaviors similar to those shown 
by the simulation. These results confirm and complete the findings of the earlier 
experimental study of the PCR’s dynamics at various Fresnel number, given in chapter 
5. The arbitrariness of choice of the modes to include in the dynamics is the main 
shortcoming of the truncated modal expansion approach. Nevertheless it appears to 
be a convenient and computationally nondemanding tool to use to study and predict 
qualitatively the spatiotemporal behavior of phase-conjugate resonators and perhaps 
of other nonlinear optical devices as well. 
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7.1 Discussion and summary 

The material presented in this thesis is the result of an inquiry about the possible 
use of a phase conjugate resonator (PGR) as an image storage as well as a processing 
element in structural control (e.g. for the real time monitoring of the deformation 
of large structures). The dynamics of the PCR has been studied in two distinct 
regimes. Below threshold (gain < loss), the PCR cannot self-oscillate. In this regime, 
the output beam is always stationary in space and time, after the transients have died 
out. This stability makes it possible to use the PCRs in image processing applications. 
When the gain threshold is exceeded (above threshold), the PCR self-oscillates and 
spatiotemporal instabilities in the output beam set in. 

This thesis was thus organized in two parts. In the first part, we have investigated 
the transient behavior and stability properties of a PCR below threshold. In the sec- 
ond part, we have studied experimentally and theoretically the PCR’s spatiotemporal 
dynamics above threshold. 

We have presented the transient dynamics of a PCR below threshold based on 
time-dependent coupled wave equations for FWM in a photorefractive crystal with 
two distinct interaction regions caused by the feedback from a slightly tilted ordi- 
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nary mirror. First, we have solved the steady state equations for the cavity s fields 
analytically, assuming nondepleted pumps and an absorption free medium. The so- 
lutions predict that the PCR self-oscillates (reaches threshold) in a broader range 
of values of pump ratios if the PCM has a higher values of the coupling parameter. 
The PCR’s transfer function (below threshold) is obtained through applying a sim- 
ple frequency- domain transformation technique to the steady state equations. The 
transfer function is shown to be largely dependent on the coupling parameter and the 
pump ratio. Contour plots of the transfer function as a function of complex angular 
frequency indicate that there is no sign of instability before the PCR can self-oscillate, 
while unstable poles quickly develop once the PCR is above threshold. 

In the transient regime of the PCR, we have successfully applied a numerical 
scheme for integrating the photorefractive FWM equations separately in the space 
and the time domains using the adiabatic elimination process. Numerical calculations 
of the cavity’s buildup and decay times as well as of the specularly reflected and 
phase conjugate reflected intensities as functions of a number of system parameters 
confirm that the parameters affecting the threshold most directly are the coupling 
parameter and the pump ratio. We have presented results which indicate that the 
cavity buildup and decay times can be tailored by varying several system parameters. 
Results predict the possible use of the PCR as an image processing element. For 
example, the PCR could, in principle, perform constrast enhancement during buildup 
and contrast reduction during decay. This same feature could be used to implement a 
time varying spatial filtering if the projected field in the PCM is the spatial spectrum 
of an image rather than the input image itself. Experimental measurements of the 
cavity’s buildup and decay times confirm qualitatively the result of the numerical 
simulations. 

In the above-threshold regime, we have studied experimentally some aspects of the 
spatiotemporal behaviors of the PCR as a function of the Bragg detuning bd achieved 
by misaligning one of the two pump beams, and as a function of the degree of trans- 
verse confinement controlled by varying the resonator’s Fresnel number F. Temporal 
aspects of the beams’ complexity were revealed by analysing the local intensity time 
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series, power spectra, and embedded phase space portraits. In the case of varying F, 
we observed bifurcation phenomena typical of low dimensional dynamical systems, i.e. 
in order of increasing control parameter, a bifurcation from a stationary to a periodic 
state, followed by more subharmonic bifurcations, the occurrence of a chaotic state, 
the occurrence of intermittency, and finally a fully turbulent state. For a particular 
case of irregular oscillations, the correlation dimension of the attractor constructed 
from the embedded phase space was calculated to be approximately 5.2, indicating 
that the attractor may be strange. Furthermore, a positive finite lower bound value 
of the Kolmogorov entropy can be extracted from the same time series, indicating the 
existence of a deterministic chaos. In the case of varying bd, the temporal aspects of 
the dynamics in the local time series range from simple periodic oscillations, followed 
by quasiperiodic states interrupted by chaotic states. All experimental results thus 
far have consistently revealed the presence of two distinct frequencies in the power 
spectra in some regions of parameter space. 

Regarding the spatial aspects of the beams’ complexity, we have demonstrated the 
existence of wavefront dislocations or vortices in the optical fields by video recording 
and by interferometry. A rich variety of defect movements has been recorded and 
analysed. The spatial complexity was characterized by the number of optical vortices 
and their variety of movements, as well as by their relationship with the beam’s 
spatial coherence. The number of defects and the spatial correlation index have been 
shown to follow similar trends as the transverse confinement of the system is relaxed, 
pointing to the role played by the wavefront dislocations in the dynamics. 

We have also tested a model based on a truncated modal expansion of the cavity 
modes, to investigate the spatiotemporal dynamics of the PCR. We have clearly shown 
that the optical vortices in the wavefront are also present in the numerical solutions. 
Simulations using the Bragg detuning as a control parameter have shown that optical 
vortices begin to appear under conditions of increased local gain which appears to be 
obtained by the small misalignment of one of the pump beams from the Bragg angle. 
Power spectra of simulated time series also contain two distinctive frequencies for a 
certain range of the off-Bragg parameter, indicating that the dynamics is governed 
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by two different time scales. The short one, on the order of the material’s damping 
time, seems to correspond to the recovery rate of the local wavefront after a vortex 
has moved through it, while the long one is associated with the alternation from one 
mode to another in the competition process and seems to correspond to the full cycle 
of the vortices’ nucleation, motion and annihilation. Numerical results also indicate 
that the loss of temporal coherence and the onset of temporal chaos in the local 
intensity fluctuations is likewise correlated to the loss of spatial confinement of the 
vortices’ trajectories and the spatial coherence. 

Thus far, several reports have been published in the literature concerning tem- 
poral instabilities in photorefractive oscillators. Each of these systems has provided 
a physical picture of how chaotic oscillations can occur. Our study shows that the 
spatiotemporal chaotic behaviors in a PCR may possibly be attributed to the mech- 
anisms of the competition among various cavity modes and the gain enhancement 
which arises from off- Bragg FWM. As we have seen in chapter 2, the steady state 
solution predicts that the range of pump ratios R p , within which the feedback param- 
eter /? can be greater than one, becomes wider and wider as the coupling parameter 
'y a d is increased. The direct numerical integration of the FWM equations in chapter 3 
confirms that the PCR can fall into the above- threshold regime in a wider range of 
pump ratios with increasing values of the coupling parameter. In chapter 6, we pre- 
sented the three dimensional FWM equations in order to account for the transverse 
dynamics in the PCR. We introduced a modal decomposition method in which cav- 
ity fields are basically described by a set of participating cavity modes. Numerical 
results show that chaotic oscillations set in for a large enough coupling parameter 
(e.g. | 7 0 d| >8). In the modal decomposition method, we have seen how the modal 
amplitudes of the cavity fields are coupled to each other in eqs. 6.20—6.23. It is this 
complex coupling that brings about the unstable oscillations which would otherwise 
not be seen in the plane-wave FWM model even in cases of quite large coupling pa- 
rameters. The off-Bragg parameter bd obtained by tilting one of the two pumps away 
from the Bragg angle is another significant factor in lowering the threshold for the 
instabilities observed in PCR. When bd = 0 (no Bragg detuning), only periodic os- 
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dilations are seen in the output. When bd ^ 0 (with Bragg detuning), quasiperiodic 
and irregular oscillations appear in the output. 

With no Bragg detuning, the coupling between cavity modes provides only chan- 
nels for energy exchange between modes. This “beating” between oscillating modes 
leads to periodic changes in the output beam. With Bragg detuning, there is an 
additional phase transfer between the interacting modes due to the additional phase 
shifts of the individual mode gratings. It is this coupling of the phases which makes 
possible the occurrence of chaotic behavior. An intuitive explanation for this behavior 
is the following. Consider what happens inside the crystal as the FWM process pro- 
gresses in time. Initially pump 2 illuminates the crystal and at time t = 0 pump 1 is 
turned on and scattered noise in the cavity forms beam 4. Then, beam 1 and beam 4 
start to write various mode gratings and beam 2 is scattered from these gratings and 
the scattered beam is fed back from the dielectric mirror. However, as the gratings 
start to form, they change the phase of the individual mode. This in turn causes a 
new mode grating to be written, which again changes the phases of the participating 
modes, and so on. This constant change in phases of the interacting beams leads to 
constantly writing new gratings and washing out old gratings. It is these dynamic 
gratings together with the feedback of the cavity that gives rise to chaotic oscillations 
once some gain threshold is exceeded. 

A final remark should be given. The arbitrariness of choice of the participating 
modes is the main shortcoming of the truncated modal expansion approach used in 
this work, despite the fact that justifications are given for truncating the series at a 
mode index of the order of the cavity Fresnel number. Yet, the method is shown to 
be a particularly convenient and computationally nondemanding tool for study of the 
spatiotemporal behavior of the PCR and perhaps of other nonlinear systems as well. 
Thus far, experimental results are supportive of the validity of this approach. 

In conclusion, this the first time to the author’s best knowledge that the spa- 
tiotemporal dynamics has been both experimentally and theoretically explored in a 
photorefractive PCR, and also that direct physical evidence of optical vortices has 
been unveiled in detail in the wavefronts of the optical fields from the PCR. Our re- 
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suits should prove useful in the study of photorefractive devices and in the assessment 
of their practical usefulness. They may also add a small contribution to the growing 
area of research in nonlinear dynamics in general. 


7.2 Possible future investigations 

One of the interesting properties of a PCR is that a great number of transverse 
modes is allowed to oscillate if the transverse extent of the cavity is opened up. This 
means that an optical processing device utilizing a PCR can have a high spatial 
resolution. Yet, the spatiotemporal instabilities due to the transverse mode com- 
petition reported in this thesis may impose limitations on the spatial and temporal 
performance of the device. To make a device pratical, methods for eliminating these 
instabilities are definitely needed. Some efforts have been recently made to alleviate 
the spatial and temporal instabilities in a PCR with BaTi0 3 as the PCM, by using 
an intracavity mode homogenizer such as a multimode fiber taper or light pipe [107]. 
The demonstrated improvement in the spatial and temporal stability, and also the 
overall increase in the transverse extent of the oscillating mode, will undoubtedly 
make possible many new processing techniques based on PCR in the above threshold 

regime. 

The mode competition process in the PCM may be favorable for other applica- 
tions. For example, associative recall and processing of images using photorefractive 
oscillators have already been demonstrated [9, 10, 11]. In these applications, compe- 
tition between modes is fundamental for the discrimination capability of the memory. 
A recent investigation to explore the possibility of using a laser as a nonlinear ele- 
ment in an associative memory architecture [108] may be an important new approach 
utilizing the continuous field distributions in a cavity rather than a discrete array of 
elements as in many neural networks. One basic requirement for this approach to 
be successful is the coexistence of many stable states, or “spatial multistability . In 
isolated multistable systems the final state is completely unpredictable because ini- 
tial fluctuations determine its evolution. Yet, different behavior is expected when the 
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system is driven externally by an identical system. In this case, a weak signal from 
the driving system is injected into the system while it is below threshold. As soon 
as the injected signal is switched off, the main system is switched above threshold. 
The system then starts with the seed left by the signal, and the evolution toward a 
stationary state is no longer random but is a function of the spatial configuration of 
the initial seed. Still, more demonstrations of practical architectures of such kinds are 
definitely needed. It seems that a PCR with high mode volume may be a candidate 
for the nonlinear element of an optical associative memory. Nevertheless, the first 
question to answer is whether spatial multistability can be achieved in systems based 
on photorefractive oscillators. To this end, the PCR patterns must be stabilized and 
simultaneously some process must be found to switch from pattern to pattern. Per- 
haps the example of using an intracavity mode homogenizer in PCR is a first step 
toward achieving these goals. 

Another aspect of the spatiotemporal dynamics concerns the behaviors of the cav- 
ity fields in other PCR geometries. The similarities and differences obtained with 
different setups may be important. At least, with such a comparison we could have 
some ideas of what is generic or geometric-specific in the observed dynamics. Fur- 
thermore, other new methods for the study of extended systems such as subspace 
Lyapunov spectra. [84] could also be applied to the diagnosis of experimental data 
from the PCR. 
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Derivation of the four- wave 
mixing equations 


In this appendix we derive the set of coupled wave equations for waves propagating 
in the PCM that are used in studying the dynamics of the PCR. At first we consider 
the optical fields as uniform plane waves of infinite extent, then we allow the optical 
fields to have transverse spatial variations. 


We start with the Maxwell’s equations for the macroscopic 
In general these equations take the form. 

variables of the fields. 

i m 

V x £ = -zr, 

C ot 

(A.l) 

-> 1 dV 47 r z* 

V X H = ~~n7 4 3 ’ 

c ot c 

(A.2) 

V • V = Anpj, 

(A. 3) 

V-B = 0, 

(A.4) 

with 


T> = £ + 4ttP = e 0 (e r0 + Ae r )^, 

(A.5) 

H = B - 4ttM, 

(A.6) 
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J = <?£, (A. 7) 

where £ is the total optical electric field, D is the electric displacement, T* is the 
induced polarization, H is the magnetic field, B is the magnetic induction, AA. is the 
magnetization, pj is the free charge density, J is the current density, c is the speed 
of light in vacuum, e 0 is the permittivity in vacuum, e TO is the background (average) 
relative dielectric constant when only the pump beams are present, Ae r is the relative 
induced dielectric constant which is slowly varying in both time and space when both 
the pump beams and the probe beams are present, and a is the conductivity. 

For ferroelectric oxides, H = B and p f = 0. Taking the curl of eq. A.l and 
combining with eq. A. 2, we obtain the wave equation: 

2( ? td 2 £ _ At d 2 £ ina d£ 

- + ( A-8 ) 

where t = t 0 t TO and Ac(r, t) = e 0 Ae r . It is assumed that Ae(r, t) is smooth in space 
and time. Suppose that the fields inside the PCM are uniform plane waves of the 
same polarization, say the extraordinary polarization (see figure A.l). We also assume 
that all waves are propagating with small angles to the z-axis and the total intensity 
profile as well as the resultant refractive index vary periodically in x with frequency 
k x . Thus, we can write down the scalar total optical field in fourier expansion as 

OO 

£{x,z,t)= J2 Mz,t)e i l ut * k ’* +,k **\ (A. 9) 

/= — OO 

and the induced dielectric constant as 


Ae(x, z, t ) 


e i( z ,t) e 


ilk x x 


(A. 10) 


with k x — k 0 y/tsmQ , k z = k oy /lcosd, k 0 = 2n/\ 0 = cu/e, where A 0 is the light 
wavelength in vacuum. 


In the following, we use the slowly varying envelope approximation (SVEA) which 
assumes that Ai(z,t) and ei(z,t) are smooth functions of 2 and t , i.e. 


d 2 A t 

dt 2 


< 2u 


dAi 

dt 


and 


d 2 A, 

dz 2 


< 2k z 


dA\ 

dz 


(A. 11) 
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ft = t 



Figure A.l: Schematic diagram for the four- wave mixing geometry. Aj and A 2 are 
the two external pumps. A\ is the probe beam and A 3 is the phase conjugate beam. 
Arrows show the beams’ polarization and direction of propagation. Beams are sym- 
metrical about Z and 6 is the angle inside the crystal, c. denotes the crystal’s axis. 
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Substituting eq. A. 9 and A. 10 with A. 11 in A.8 and equating the /th component of 
the optical field on both sides, we obtain 


dAi we dAi i(l — l 2 )kl 2n ua 

~dl cNT z 1h 2 ~k z Al± ~^kT Al 


ik 2 


(A. 12) 


In four-wave mixing, without Bragg missmatch, beam 1 is collinear with beam 2 and 
beam 4 is collinear with beam 3 (see figure A.l), we have 


k\ k 2 — &3 -f- k 4 


(A.13) 


The total optical field as a sum of four beams takes the form 

£{x,z,t) = J2A J (z,t)e^ t - n ^ (A. 14) 

j = 1 

Next, for the Bragg’s condition in a thick volume hologram situation, i.e. A; = 0 for 
|/| > 1 and A±i ^ 0, eq. A. 12 becomes 


dA± i ( we dA± i _ ik 2 0 A a A 

dz + c 2 k z dt - 2 k z e±2ATl COS0 A±U 

for waves propagating to the right, and 

dA ± i u>t dA±i _ ( ik 2 a A ( a A 

dz c*k, dt _+ 2 *, £±2 fl + cos 


(A. 15) 
(A-16) 


for waves propagating to the left, where a = ‘l-Kajcyft. Thus we can write down the 
coupled differential equations for waves propagating inside the PCM in the notation 
of figure A.l as follows: 


dAi 1 dAi —ik 0 a 

dz ^ v cos 0 dt 2 v /ecos0 C2 4 cost? 1? 

dA 2 1 dA 2 ik Q a 

dz v cos 0 dt 2 v /ecos0 2 3 + cos0 2 ’ 

dAz _ 1 dA 3 _ ik 0 a 

dz vcos0 dt 2^/tcosO 2 2_ ^cos0 3 ’ 

dA 4 1 dA 4 _ -ik a a 

dz A vcos0 dt 2 v /ecos0 C_2 1 cos0 4 ’ 


(A.l 7) 

(A.18) 
(A. 19) 
(A. 20) 
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with v = cjy/t. From eq. A. 14, the total intensity is 


I{x,z,t) = \£\ 2 

= I 0 (z,t) + {A l Ale~ i ^ )? + A 2 AZe- i &- k) - f 

A;A 4 e-^- ?2) 0 + c.c. 

(A. 21 ) 

with I 0 {z,t) = Ej=i |Aj(2,<)| 2 . Because of the small incident angle to the z-axis for 
all fields, the exponentials in the 4th and 5th terms are fast changing phase factors 
with respect to the slowly varying amplitudes so that they becomes insignificant in 
space averaging. Now we let k\ — = k 2 — = Kx with \K X \ ~ 2 k x . Then the light 

intensity can be written as 

I(x,z,t) = I 0 {z,t) + [h{z,t)e~ iK * x + C.C.] , 

(A. 22) 

with 

I 2 i_z.fi) = Ai[z, t)A 4 (z : t) 4- A 2 (z, t) A^iz ,£) . 

(A.23) 

The induced dielectric constant is related to a space-charge field via 
electro-optic effect (see Appendix B) and is given by 

the linear 

Ae(^, 2 ,f) = — ■\/th 0 r e j jE a dx , z, t), 

(A.24) 

where n 0 is the bulk refractive index, r e // is the electrooptic coefficient, and E sc is 
the space-charge field determined by the transport equations B.1-B.4. Also we can 
expand the space-charge field in the same form as Ac in eq. A. 10: 

E ac = E 0 +[E 2 {z,t)e- iKxX + c. c : . 

(A. 25) 

Thus, comparing eq. A. 25 with A. 10 and equating the coefficients of e 
sides of eq. A. 24, we have 

xX on both 

e +2 (z,t) = ^Yh 3 0 r eff E 2 (z,t ), 

(A. 26) 

and 

e-a(M) = c + 2 ( 2 >*)- 

(A. 27) 
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Taking 9 < 1, and using eqs. A.17-A.20 with eqs. A.26 and A. 27, 
of coupled wave equations having the form: 

we obtain a set 

dAx 1 dA\ 

dz + v dt = -QA,-aA u 

(A. 28) 

dA 3 1 dA X 

dz -- v dt -- QA ^ aA - 

(A. 29) 

9A 3 1 dA 3 „ , . , 

dz v dt - QA2 + qA3 ' 

(A. 30) 

dA\ 1 dA * „ 

17 + v dt 

(A.31) 

with the index grating 

„ icon 3 „ 

Q — ^ c r effE 2, 

(A. 32) 


where E 2 is given in eq. B.32. 


Next, we consider the possible transverse spatial variations of the optical fields. 
It means that the fields’ amplitudes depend in general on x and y as well as z, i.e., 

4 

s (x,y,z,t) = '52 A j(x,y,z t t)e* ut - S rt + c.c. (A.33) 

3=1 

The wave equation A. 8 can be written as 


V 2 £(x,y,z,t) - 


(A. 34) 


where n 2 - e + Ae. Also when considering weak modulation in a photorefractive 
medium, we can expand n in the same form as the induced dielectric constant: 


n (x,y,z,t ) = n(x,y,z) + n 1 e'^m(x,y, z, t)e'^ r + c.c. 

= h + An(x,y,z,t), 

where m(x,y,z,t ) is the first order perturbation which is given by the ratio 

m<», »,».*) = 

I 0 (x,y,z,t) 

_ A1A4 + A^As 
To 
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and n is the bulk refractive index, An is the induced refractive index. Here I< x = 
k 4 - fa = fa - fa. We assume a small modulation in the refractive index, i.e., 
| An | <C h. Then the wave equation A. 34 can be reduced to 


h 2 


V 2 £(x,y,z,t)- — + 2n 


An\ d 2 £ toadS 


dt 2 


dt 


= 0. 


(A.37) 


Substituting eqs. A. 33 and A. 35 in eq. A.37 and using the SVEA, we obtain 


[v£ - kf - 2 iki • V 


Ai - 


i=i 


- 2 




-me 


+ 2 "4) A 


~ ( J) + A AJe-* '' + c.c. = 0, ( A.38) 


where 


2 _ d 2 ^ 

Vy = tt - - + 


, (A. 39) 

dx 2 ' dy 2 

is the transverse Laplacian operator. Equating the coefficients of e '** r separately, we 
obtain the following set of equations for waves propagating inside the PCM including 
the transverse spatial variations in the fields’ amplitudes: 


'»V 


d 


Ai — —QA 4 , 


pi +k '- v+a+ M, 

' iX7 2 

T + fa • V + a + 


2k 


vdt 


a 3 = -qa 3 . 


■w +i4 ' v+a+ J<) 45 


QA\ 


(A.40) 

(A.41) 

(A.42) 

(A.43) 


where we have ignored the products of two small quantities viz n\ dAi/dt as well as 
the terms (1 /u)dA,/dt; a = 2t ra/c^/l and v = c/y/i. The grating Q is given by 


Q{x,y,z,t) = 7 0 m(x, y, z, t), 
where the coupling coefficient -yo I s 


7o = 


toon ie 


,i<t> 


2c 


(A. 44) 


(A. 45) 


The expressions for the amplitude m and the phase of the refractive index modu- 
lation are given in eqs. B.55 and B.56 of Appendix B. 
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The Kukhtarev equations 


In this appendix we use Kukhtarev band transport model [42] to derive the ba- 
sic equations governing the hologram formation process in photorefractive media. In 
Kukhtarev model, the photorefractive materials are assumed to contain donor and 
acceptor traps which arise from the crystals’ impurities. These traps introduce inter- 
mediate electronic energy states in the bandgap of the insulators. When photons of 
sufficient energy are present, electronic transitions by photoexcitation take place. As 
a result of the transitions, charge carriers are excited into the conduction band and 
the ionized donors become empty trap sites. Furthermore, movements of the charge 
carriers are generally described by diffusion, drift and the photovoltaic effect [16]. 
The migration of the charge carriers from bright regions where they are photoexcited 
to dark regions where they are recombined with acceptors builds up a space-charge 
field inside the crystal. The space-charge field modulates the material’s refractive 
index via the linear electrooptic effect [40] (Pockel effect). 

The space-charge field created by the migration and trapping of the charge carriers 
is determined by the following set of transport equations: 

dn dNt 1 

Yt = ~ q V ' 3 (continuity), (B.l) 

— = (SI + (3)(Nd — N£) — 'y R nN£ (rate equation), (B.2) 
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j = qpnE - k B TpVn + ple c (current equation), (B.3) 

V • E = - (n + Na - (Poisson’s law), (B.4) 

where (SI + / 3)(N D — Np) is the rate of charge carrier generation, 7 R nNp is the rate of 
trap capture. Here S is the cross-section of photoionization, /? is the rate of thermal 
generation, 7 * is the carrier-ionized trap recombination rate coefficient and n and 
Np are the concentrations of the carriers and ionized traps, No is the concentration 
of the donor traps, I is the light intensity, N A is the acceptor concentration, p is 
the mobility, T is the temperature, k B is the Boltzman constant, e is the dielectric 
constant, pi is the photovoltaic current and p is the photovoltaic constant, e c is the 
unit vector along the c axis of the crystal, E is the total field which includes the 
induced space-charge field E sc and any external or internal fields such as internal 
ferroelectric fields; and q is the charge of the charge carrier. 

In general p and e are tensors and N D and N A could vary with position. For 
simplicity we will assume that p and e are scalars and N B and N A are constant 
throughout the material. Eqs. B.1-B.4 are nonlinearly coupled and generally difficult 
to solve without some approximations. In the following we will simplify the band 
transport equations using various approximations to solve for the space-charge field 
and obtain the basic equation that governs it. 

We begin with a perturbation scheme to linearize the above coupled equations. 


To this end, we let (see figure B.l) 

I(x,z,t) = I 0 (z,t) + \h(z,t)e- iKlX + c.c.] (B.5) 

1 F(x, 2 , t ) = E a (z, t) + [fi 2 ( z > t)e~ ,KxX + c.c.] , (B.6) 

n(x, z , t ) = n 0 (z, t ) + [« 2 (^i t)e~ lKxX + c.c.] , (B.7) 

Np(x,z,t) = Np 0 (z, t) + [N£ 2 (z, t)e~' hxX + c.c.] , (B.8) 

j(x,z,t) = jo(z,t)+ [j 2 (z,t)e-' KxX + c.c] . (B.9) 


Here E(x , 2 , t ) is the overall electric field in the photorefractive material, which is the 
sum of any external field E 0 and the induced space-charge field E 2 , and we assume 
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Jt = d 





Maximum of the light intensity 
Maximum of the grating 


Figure B.l: Hologram formation in a photorefractive material using the four-wave 
mixing scheme 
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that 1 2 ( 2,0 varies relatively slowly along the z-axis as compared to the fast change 
of e~ ,KxX along the x-axis, etc. (SVEA). 

At first we consider the quasi-steady static case, i.e. dn Q /dt ~ dN£ a /dt ~ 0. We 
also assume that there are a lot of available donors and acceptors to linearize the 
generation and recombination rates, i,e., No No — N^ and Na > n. From now 
on, we also ignore the small thermal factor, i.e. /3 < SI. Then, from eq. B.4 and 
assuming the SVEA, we have 


^ ^ q -(n + N a - JVj). (B.10) 

ox e 


We obtain the following equations upon equating the zeroth order terms and the 
coefficients of e~' hxX : 



«D„ = + n, ~ N Al 

feu) 


Ei = #-(n 2 - N&). 

(B.12) 

Consider a maximum space-charge field with charge separation comparable to the 

grating period such that 

K x E' q = 

(B-13) 

Then, in eq. B.12 

iE' 

E 2 = n2 - AT+ 2 ). 

rsD 

(B.14) 

Now, from eq. B.2, we have 



9Nd 2 „-iK x x _j_ n n ^ 

dt + 

S [l 0 + I 2 e~ iKxX + C.C.] N d 

(B.15) 


- 7 r [«o + n 2 e~ thxX + c.c.] 

[■'Vd. + + c.c.] , 

which gives us 

SI 0 N d SI 0 N d 
n ° " IrNdo ' ’ 

(B.16) 
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For the quasi-steady state, or when the generation rate equals the recombination rate, 
i.e. dN^ 2 /dt ~ 0, we have 


/v+ 

iV L>2 


S 1 2 Ad N a 

n 2 . 

lRn 0 n 0 


(B.18) 


Combinding eqs. B.l and B.3, we get 
I" dn 


dt 


2 e~ iKxX + c.c. 


-iKrX , 

~dT e + c - c - 

-fi In 


(B.19) 


■P 

k B Tfi 


o + n 2 e ,AlX -f c.c.J [(-i7C)£ 2 e tAlX -f c.c.j 
(-*7C)n 2 e- ,A '^ + c.c.] [E 0 + E 2 e~ iK * x + c.c. 


[(— z7l x ) 2 n 2 e _,AlX + c.c.j 
a L J 


_P 

As dn 2 /dt — dN^/dt ~ 0, we get 

(e 0+ - ,,TI<X 


\—iK x )I 2 e~ ihxX + c.c.] . 


«2 = - 

£2 + — (£) 

/ 

\/ o y_ 




(B.20) 


We recognize the diffusion field E D = k g TK x /q and the photovoltaic field E p = 
pl 0 /qfin 0 . Then we rewrite n 2 as 

£2 + E p I 2 j I 0 


Tin 


E 0 4- iEd n °' < B - 2I > 

Now, substituting eqs. B.18 and B.21 in eq. B.14, we have the following static space- 
charge field: 

E% = -iE> f + !l\ . (B.22) 


N 


D 


H~ Ep^/Ip I2 

E 0 + iE[j I Q 


Since 


E 


7 ( NA+n Q \ qN A q(N D - N A ) at _ qN e _ „ 

uc/vo ~ Na = 7 k : = e ” 


N, 


D 


tK x tK x N D l ' A ~eK x -" q ' ^ B ‘ 23) 

where N e is the effective trap density, we can simplify eq. B.22 such that the modu- 
lation of the static space-charge field is 


£2 = -»£, 


(£0 + E p ) + iEp 

.£ 0 + ?(£, + £: D ) J 


©■ 


(B.24) 
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Next, we examine the time development of the space-charge field. In the short 
writing time limit, i.e. shortly after illumination, N£ will not change greatly from its 
dark values (i.e. N% ~ N A )• Then eq. B.2 may be written as 

= SI(N D - N a ) - (B.25) 

ut T r 


with t r = \h R N A . Substituting this in eq. B.l, 

^ = S1(N D -N A )-f 


dt 


d ( p k B Tpdn p > 

~Tx aJ + ^. 


(B.26) 


As we expect dn/dt — 0(n/T d ;), whereas is the dielectric relaxation time, for 
illluminating times greater than a few t r we can then make the approximation that 
T R dn / dt 0. Also, for times short compared to r di , the space-charge field will not 

change greatly and we may make the further approximation E sc ~ E 0 . Then eq. B.26 
will reduce to 


t r S h + I 2 e~ iKzX + c.c.] ( N d — N a ) — n — pr R E t 


dn 

dx 


(B.27) 


k B Tp d 2 n p 

t r 7TT ” Tr ~ 

q dx 1 q 


+ c.C. 

OX 


which, after equating the coefficients of e xKxX , gives 


1 + iEjE}± n ( h\ 

1 + E d / Dm — iE 0 / Em °\I 0 ) 

_ Em + iEp n ( J2. \ 
\E m + E d ) - iE 0 \ n ° \lj ’ 


(B.28) 


with Em = 1 /{pr R K x ). In eq. B.24, as ~ JV*, then we have for short time 


n£ ~ n a + (sin d - t. 

(B.29) 

So, from eq. B.4, and with N A , N D - N A > n , we have 


dE ac = 9 | r IL_ 1 s/iv D > ) t, 
dx e \r n / 

(B.30) 

which gives, by equating the coefficients of e~' KxX , 


E 2 = n ( n2 5/ 2 A d )<. 

eAj Vt* / 

(B.31) 
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Therefore, combining eq. B.28 and eq. B.31 we have 


E 2 = 


Ed ~ i(E 0 + E p ) 


c \_E a + i(E\f + Ed) 
Em I" Ep — i(Eo + E p ) 


IoJ 


(B.32) 


r dt [E 0 + i(E\f + Ed)\ 




with T d , — e/(qfin 0 ). If we assume that the space-charge field builds up exponentially 
such that 


E 2 = E*{ 1 - e~ vt ) 

~ E^rjt 


(B.33) 


for the short time limit, then, comparing this with eq. B.31 and using eq. B.24, we 
have 


E 


M 


E q + i(E q + Ep) 


1 


1 

r d , 

1 


E 0 -f i(E\j + Ep)\ T d , 

( 1 + §)(' + fc )+ E ~ 


(e^Y 

\Em) 


EqE M 


i (E^ _ E, 

\E m E, 


) 


+ ( 1 + fe) (fc) +(> + fc)’J 


— (VR + iVi)- 


Therefore, we have a complex time constant for the material. Thus, 


e -nt _ e -t/r, e iu,t 


with 


and 


r dt _ 

i 

(! + 

MjA 

E m ) 

2 / 
1 +( 

eA 

Em) 

2 

1 

7~dt 

r 

i 

(>+t) 

1 1 

(• + & 

)+ 

E n 

EqEM 


_ Vi 




Eg Eg 

Em E q 


Ti> L(fe) + ( 1 + E^) J 

In the following we define various time scales of the material: 


T n = 


nK x E 


D 


(diffusion time), 


(B.34) 


(B.35) 


(B.36) 


(B.37) 


(B.38) 
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1 


T P = 


jX I\. x;Eq 

1 


Ti = 


SI 0 

Tr = 


1 


(drift time due to external field), 
(inverse of photoexcitation rate), 
(recombination time), 


Tj. = 


7 r n a 

1 (dielectric relaxation time). 


qp.n 0 

From the above definitions, we also notice that 

Ed 


Em 

A 

Em 

Ed 

E 0 

Em 

X 


Tr 

i 

Td 

t r 
” ; 
T E 

Te 

Td 

T* 

Ti 


(B.39) 

(B.40) 

(B.41) 

(B.42) 

(B.43) 

(B.44) 

(B.45) 

(B.46) 


Then the material response time and the grating natural frequency can be expressed 
in terms of these various time scales as 


Te 7~di 


1 

(1 + S) 

2 

1 + 

(iR 

\TE 

) 2 

(1 + (> + 

td ) 

+( 

\ 2 

ir) 1 

T E/ 1 

(7) 


and 




1 1 

\T E ) 

l( 

I0L — 

v T/ 

0 


Xt j 

fl + za 

t TD, 

. 2 

) +< 

(3 

r 


(B.47) 


(B.48) 


Here we can see that u g = 0 when E 0 = 0, i.e. the material time constant becomes 
real when there is no applied field. Now we consider eq. B.33 as a solution for E 2 of 
a differential equation of the form 


9E 2 n _ nS 
“aj b 1E2 — f]E 2 . 

at 


With eq. B.34, the above equation can be rewritten as 


dE 2 ^ 

TV,.— + 


dt 


(E m \ 

E 0 + i(E q + E d ) 

UJ 

_E 0 + i(Em + Eq) 

(E m \ 

E 0 + i(E q + E d ) 

UJ 

_E q + i(Em + Ed)_ 


(B.49) 


(B.50) 


E%. 
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Using eq. B.46 we have 


Tl 


dIX 

dt 


+ 


E a + i(E q + Ep) 
E 0 + i(Em + Ed) 
E 0 + i(E q + Ed) 
E 0 T i(E\f T Ed) 


E 2 



(B.51) 


From eq. A. 32 of Appendix A we know that the grating is directly proportional to 
E 2 . Thus we can write down the time equation for the grating as 


dQ 

dt 


+ CQ — — C(A\A * 4 + A2A3), 

n 


(B.52) 


with 


and 


where 


and 


E 0 + i(E q + Eq) 
E 0 + i{E\f + Ed) 


7 0 


4 c 
iun 1 


r eff (-iEg) 
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i<t> 


( E 0 + E p ) + iEd 
E 0 + i(E q + Ed) 




1 

~2 r ' ff 


nlE„ 


(E 0 + E P ) 2 + E 2 D 

E 2 0 + (E q + E D y 





(j> = tan 


-1 


EpjEg + Ed) + K(Eo + Ep) 

E 0 E q + E p (E q + Ed) 


(B.53) 


(B.54) 


(B.55) 

(B.56) 


When E 0 = E p = 0, <j> = w/2. Thus, is real. From eq. B.52, this means that 
when there is no external applied field, the grating does not have a natural oscillation 
frequency and it is always 7 t/ 2 phase shifted with respect to the interference fringes 
of the optical fields. In addition, the response time is just inversely proportional to 
the intensity. 
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